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In the present tutorial we address a problem with a long history, which remains of great interest 
to date due to its many important applications: It concerns the existence and stability of periodic 
and quasiperiodic orbits in iV - degree of freedom Hamiltonian systems and their connection 
with discrete symmetries. Of primary importance in our study is what we call nonlinear normal 
modes (NNMs), i.e periodic solutions which represent continuations of the system's linear normal 
modes in the nonlinear regime. We examine questions concerning the existence of such solutions 
and discuss different methods for constructing them and studying their stability under fixed and 
periodic boundary conditions. 

In the periodic case, we find it particularly useful to approach the problem through the discrete 
symmetries of many models, employing group theoretical concepts to identify a special type of 
NNMs which we call one-dimensional "bushes" . We then describe how to use linear combinations 
of s > 2 such NNMs to construct s-dimensional bushes of quasiperiodic orbits, for a wide variety 
of Hamiltonian systems including particle chains, a square molecule and octahedral crystals in 
1,2 and 3 dimensions. Next, we exploit the symmetries of the linearized equations of motion 
about these bushes to demonstrate how they may be simplified to study the destabilization of 
these orbits, as a result of their interaction with NNMs not belonging to the same bush. Applying 
this theory to the famous Fermi Pasta Ulam (FPU) chain, we review a number of interesting 
results concerning the stability of NNMs and higher-dimensional bushes, which have appeared 
in the recent literature. 

We then turn to a newly developed approach to the analytical and numerical construction of 
quasiperiodic orbits, which does not depend on the symmetries or boundary conditions of our 
system. Using this approach, we demonstrate that the well-known "paradox" of FPU recurrences 
may in fact be explained in terms of the exponential localization of the energies Eq of NNM's 
being excited at the low part of the frequency spectrum, i.e. g = 1, 2, 3, .... These results indicate 
that it is the stability of these low-dimensional compact manifolds called g-tori, that is related 
to the persistence or FPU recurrences at low energies. Finally, we discuss a novel approach to 
the stability of orbits of conservative systems, expressed by a spectrum of indices called GALI^, 
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k — 2, ..., 2iV, by means of which one can determine accurately and efficiently the destabilization 
of q tori, leading, after very long times, to the breakdown of recurrences and, ultimately, to the 
equipartition of energy, at high enough values of the total energy E. 

Keywords: discrete symmetries; Hamiltonian systems; nonlinear normal modes; periodic and 
quasiperiodic orbits; stability; chaos 



1. Introduction 

The dynamics of Hamiltonian systems, or, more generally, conservative mechanical systems (preserving 
phase space volume) is of great importance for the understanding of many problems arising in classical 
mechanics, astronomy, solid state physics, plasma physics and nonlinear optics. The simplest states of 
such systems are their equilibrium points, where all variables are fixed for all time, and their periodic 
orbits, where all particles oscillate and return to their starting values after a time interval T. The first 
question concerning these states is their existence. Once this has been established, what one needs to 
know is their stability^ i.e. their behavior under small perturbations of the initial conditions and parameter 
values. Two great contemporary scientists were the first to deal with these issues in a systematic and 
comprehensive way: The Russian mathematician Alexander Mikhailovich Lyapunov (1857 - 1918) and 
the French mathematical physicist Henri Poincare (1854 - 1917). Lyapunov devoted a great deal of his 
efforts to 'local' stability analysis, obtaining specific conditions for the beha vior of solu t ions of systems 
of ordinary differential equations (ODEs) in the vicinity of equilibrium points [Lyapunovl . 118921 ]. Poincare 



was more concerned with 'global' properties of the dynamics, like non - integrability and t he occurrence o f 
irregular (or chaotic) solutions, wandering over large domains of the available state space Poincard . Il892l |. 



Remarkably, both of them worked on periodic solutions of A^ - degree - of - freedom Hamiltonian systems 
and were fascinated by the problem of the stability of the solar system. 

In the present paper, our purpose is to review certain recent results concerning the dynamics of 
Hamiltonian systems of N degrees of freedom, whose equations of motion have the form 

^ = ^^ ^ = _^^ k = l,2,...N (1) 

dt dp^. ' dt dq^. ' 

In many cases the Hamiltonian is expanded in power series in the variables of positions and momenta, as 
a sum of homogeneous polynomials Hm of degree m > 2 

H = H2 + H3 + ... + Hm = E , Hm = Hmiqi,-.-,qN,Pi,---,PN) (2) 

{E is the total constant energy), so that the origin, qk = Pk = 0, k = 1,2, .., A^ is an equilibrium point 
of the system. However, there also many examples, e.g. possessing discrete symmetries, for which no such 
decomposition is necessary as one can work directly with the original Hamiltonian. 

Assume now that the linear equations resulting from ([1]) and ([2]) with Hm = 0, for all tti > 2, yield 
a matrix, whose eigenvalues all occur in conjugate imaginary pairs, HXk and provide the frequencies of 
the so-called normal mode oscillations of the linearized system. According to Lyapunov, if none of the 
ratios of these eigenvalues, Xj/Xk, is an integer, for any j,k = 1,2., N, the linear normal modes continue 
to exist as periodic solutions of the nonlinear system ([T]), when higher order terms Hs,!!^,... etc. are 
taken into account in ([2]). These solutions have frequencies close to those of the linear modes and are 
examples of what we call nonlinear normal modes (NNMs), where all variables {qk,Pk) oscillate with the 
same frequency (A^ =)ujj = 2tt/Tj, returning to the same values after a single maximum (and minimum) 
in their time evolution over one period Tj. 

What is the importance of these periodic solutions? Once we have determined that they exist, what 
can we say about their stability properties under small perturbations? How do these change when we vary 
the total energy E m ^? Does their loss of stability affect only their immediate vicinity or can they also 
influence the dynamics of the system as a whole? These are the questions we shall try to answer in this 
paper. 
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Let us consider, as a specific example, a mechanical system in one dimension described by the the 
A^-particle Hamiltonian 

1^1^ a ^ ^ 

^ = 2^pI+2 Z^^^'^+i ~ ^'^^^ + 3 IZ^^'^+i ~ ^^^^ + 4 ^^^k+i - Xkf = E (3) 

fc=l fc=0 fc=0 k=0 

called the FPU chain in honor of F ermi, Pasta and Ulam, who were the first to study it numerically in 
the early 1950's [Fermi et all Il955l |. They discovered certain very interesting and surprising phenomena 



concerning its dynamics, which we will discuss later in this review. For the time being, let us note that 
when this chain is studied under fixed boundary conditions (f.b.c), i.e. xq = xn+i = 0, there are many 
A^'s for which the linear normal mode frequencies satisfy the incommensurability condition of Lyapunov's 
theorem and hence they may be rigorously continued to the nonlinear regime of a 7^ and/or /3 7^ 0. 

By contrast, when periodic boundary conditions (p.b.c.) are imposed, i.e. Xj = Xj+j\[,Pj = PN+j,j = 
1,2,...,N, the linear mode spectrum becomes degenerate for all N and Lyapunov's theorem cannot be 
invoked. What do we do in that case? How do we study the existence and stability of such nonlinear 
normal modes, or NNMs? As we explain at length in sections H and HI, this is an important case where 
the identification and analysis of the system's discrete symmetries turns out to be of great relevance to the 
understanding of the dynamics of the problem. 

Thus, we start by demonstrating in section H precisely how one can use the powerful techniques of 
group theory to establish the existence of families of such NNMs for a variety of mechanical systems, in- 
cluding particle chains in one dimension (under p.b.c), as well as certain 2-dimensional and 3-dimensional 
structures. Then, in section HI, we show how one can combine such periodic solutions to form "bushes" of 
quasiperiodic orbits and exploit the symmetries of the system to simplify the variational equations about 
these orbits and study the motion in their vicinity. Ultimately, of course, one would like to be able to obtain 
invariant manifolds on which the dynamics of the system is as simple as possible. This is not easy to do 
in general, but if the equations of motion possess discrete symmetry groups, one can single out some of 
these manifolds using r egular group theoretical methods developed in Sakhnenko Sz Chechinl . Il993l . Il994l : 



Chechin Sz Sakhnenkd . 119981 ]. These methods provide, in fact, the mathematical basis of the theory of 



bushes of NNMs, which plays an important role in many of the results discussed in this review. 

Before discussing this theory, however, let us illustrate its main ingredients on a seemingly simple 
example of an A^ degree of freedom Hamiltonian system. In particular, let us consider the famous Fermi- 
Pasta-Ulam /3-chain (FPU-/3), i.e. Eq. (|99p with a = 0, representing a one-dimensional lattice of unit masses 
coupled to each other by identical nonlinear springs (see Fig. [T]). The ordinary differential equations (ODEs) 

Fig. 1. FPU chain. 

describing longitudinal vibrations of the FPU-/3 chain can be written in the form 

Xi = f{xi+i-Xi) - f{xi-Xi-i), i = l..N, (4) 

where Xi{t) is the displacement of the ith particle from its equilibrium state at time t, while the force 
/(Ax) depends on the spring deformation Ax as 

/(Ax) = Ax + /3(Ax)3 (5) 

For simplicity, we assume that /3 > 0, since in that case the potential is positive definite everywhere and 
particles cannot escape to infinity. If we also impose p.b.c. we must require: 

XN+lit) ^ Xl{t), Xo{t) = XNit). (6) 
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Thus, we shall study possible dynamical regimes of this chain by attempting to solve (jl]) for the 
"configuration" vector 

X{t) = {xi{t),X2{t),...,XN{t)] (7) 

whose components are the individual particle displacements. 

As is well-known, it is pointless to try to obtain this vector as a general solution of ([3]). We will, there- 
fore, concentrate on studying special solutions represented by the NNMs mentioned above. In particular, 
let us begin with the following simple periodic solution 

X{t) = {x{t), -x{t),x{t), -x(t), ...,x{t), -x(t)} (8) 

which is easily seen to exist in the FPU chain with an even number of particles {N mod 2 = 0). This 
solution is fully determined by only one arbitrary function x{t) and is called "vr-mode" or boundary zone 
mode. 

From Eq. ([8]), it is evident that every two neighboring particles possess equal in value but opposite in 
sign displacements at any instant t. To investigate the dynamical regime of this mode in the FRU-/3 chain, 
let us substitute X{t) from Eq. ^ in Eq. ([H), taking into account Eqs. ([5]) and Q. As a result, all N 
equations ^ become identical and can be written in the form: 

,2. , ,Y/^^.3 



x + u;'x + jl^j^jx-' = 0, (9) 

where cu^ = 4, 7 = 16 and the Xi have been scaled by a factor \/yjN . 

This is the well-known Buffing's equation. In the case /3 > it describes periodic vibrations with 
arbitrary fixed amplitude A determined by initial conditions which have the form 

x(0) = A, x(0) = 0. (10) 

As is well-known, the ana lytical solution to Eqs. ([9]), (jlOp can be expessed via the Jacobi elliptic cosine 
Abramowitz fc Stegunl . ll965l |: 

x{t) = Acn{^t,k^), (11) 

where 

n'^ = u'^/{l-2k^), (12) 

while the modulus k of this elliptic function is determined by the relations 



a;2 + 6A2' '\N. 

Thus, the so-called vr-mode defined by Eq. ([8]) describes an exact dynamical regime of the nonlin- 
ear system. It re presents an example of the type of nonlinear normal modes introduced by Rosenberg 



Rozenbergl . Il962t |. From Eqs. (IIH [T3]l it can be seen that in the low-amplitude limit {A — )• 0) the above 
NNM is an exact continuation of the linear normal mode (A; — t- 0, fi — )• w), as the elliptic cosine tends to 
the ordinary cosine. 

NNMs exist in Hamiltonian systems with rather specific interparticle interactions, as, for example, in 
cases where the potential energy is a homogeneous function of all its arguments. In many cases, however, 
the existence of NNMs may be established by certain symmetry related arguments, whence we refer to 
such dynamical objects as symmetry- determined NNMs. 

The vr-mode ([8]) clearly represents a dynamical state of the form: 

X(t) = x(t){l,-l,|l,-l,|...,|l,-l}. (14) 

Hence, the following question naturally arises: Are there any other such exact NNM s in the FPU-/3 chain , 



unde r p.b.c . ? Some example s of such modes a r e already known i n the literature iPoggi fc Ruffo 



1997 



Rinkl. I2OO3I: IShinoharal . I200I l2003l : lYoshimural . l2004l : iLeo fc Leo I . l2007l : lAntonopoulos fc Bountii . 120061 
Kosevichl . 1 19931 ] . under a terminology that differs among authors. Below we list these exact states in detail 



^ 27 

X(t) = a;(t){l, 0,-1,11,0,-1, 1,..., |1, 0,-1}, w^ = 3, 7 = ^ (A^ mod 3 = 0). (15) 
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27 
X(t) = a;(t){l, -2,1, |1, -2,1,1,.. .,|1, -2,1}, u^ = 3, 7=y {N mod 3 = 0). 



(16) 



X(t) = 2;(t){0,l, 0,-1,10,1,0,-1, 1,..., |0,1, 0,-1}, a;2^2, 7 = 4 (iV mod 4 = 0). (17) 



X(t)=x(t){l,l,-l,-l|l,l,-l,-l|,..., 11,1,-1,-1}, 
a;2 = 2, 7 = 8 (TV mod 4 = 0). 



(18) 



X(t) = x(t){0, 1, 1,0, -1,-1|0, 1,1,0,-1,-11,. ..,10,1,1,0, -1,-1}, 



(19) 



1, 



I {N mod 6 



a;" = i, 7 = ^ (,-/V fTiod = 0). 

Let us comment on certain properties of the above NNMs for the FPU-/? chain, having an appropriate 
number of particles in every case: 

PI) Displacement patterns of NNMs ()14p - (J19p are defined in terms of the vibrational state primitive 
cells (divided by vertical lines), whose number of elements ttt. is called the multiplicity number. 

For (|14p - (|19p these numbers have the values: m = 2 for Eq. p4p . m = 3 for Eqs. (|15p and (|16p . 
?7i = 4 for Eqs. (fT7|l and (fT8]l and m = 6, for Eq. (fT9]) . Note that, due to p.b.c, the FPU-/? chain in its 
equilibrium state is invariant under translation by the interparticle distance a, while in its vibrational states 
it is invariant under a translation by ma. Moreover, NNMs (I14p - ()19p do not differ only by translational 
symmetry (m,a), but also by some additional symmetry transformations which are discussed later in this 
paper. 

P2) Each of the NNMs ()14p - (I19p depends on only one function x{t) and, therefore, describes a single 
one- dimensional dynamical regime. This function satisfies Buffing's equation ([9]), with w^ and 7 that differ 
according to the displacement pattern of each NNM. 

P3) Every NNM (llSp - (I19p possesses some "dynamical domain" of its own. For example, consider the 
mode (jisp : It is easy to check that cyclic permutations of each primitive cell in its displacement pattern 
produces other modes: x(t){l,0, — 1, |1,0, — 1, |..., |1,0, — 1} — >• x(t){0, — 1, 1, |0, — 1, 1, |..., |0, — 1, 1, } — >• 
x(t){ — 1, 1, 0, 1 — 1, 1, 0, |..., I — 1, 1,0, }, which differ from each other only by the position of their stationary 
particles. 

As a consequence, all dynamical properties of these NNMs turn out to be equivalent. In particular, 
as we shall show in detail in section III, they possess the same stability regions. Because of this reason, 
we only need to study one representative of the set of equivalent NNMs. Such sets are called "dynamical 
domains", borrowing the term from the theory of phase transitions in crystals. 

Many aspects of exis t ence and stability of NNMs have been discussed in the l i terature (see, for exam- 
ple. IChechin et a/.l. l20oi. l2005l: iRinkl. 120031: IPogd fc RuflFol.ll997l:IShinoh^.l200ll2003l:lYoshimu"ral .| 2004 



Leo k Leo I. l2007l : lAntonopoulos fc Bountisl . I2OO6I : iBudinskv fc Bountisl . Il983l : ISanduskv k Paed . Il994l : 
Kosevichl . 119931 ]. What is important at this point is to pose certain fundamental questions concerning 
such NNMs, which we shall proceed to answer in this review to the best of our ability: 

Ql) Is the hst (HH) - dill) of NNMs for the FPU-/3 chain complete? 

Indeed, at first sight, it seems that many other NNMs can exist, for example, modes whose multiplicity 
number m is different from those listed above. 

Q2) What kind of NNMs exist in nonlinear chains with different interactions than those of the FPU-/? 
chain? 

In fact, in most papers (see Poggi fc Ruffol . 119971 : IShinoharal . 12002, 120031 : lYoshimural . l2004l : iKosevichl . 

19931 1) the NNMs listed as (fH|) - (fT9]) above have been discussed by analyzing dynamical equations which 
are only connected with the FPU-/? interparticle interaction. 

Q3) Do there exist NNMs for Hamiltonian systems which are more complicated than the monoatomic 
chains? For example, one can pose this question for diatomic nonlinear chains (with particles having 
alternating masses), 2-dimensional (2D) lattices, or 3D crystal structures. 

Q4) Finally, there is a more subtle issue: How can one construct exact multidimensional dynamical 
regimes in nonlinear A^-particle Hamiltonian systems and study their stability by locating the unstable 
manifolds of the corresponding quasiperiodic orbits? 
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Note that the NNMs discussed so far represent exact one- dimensional regimes, because they are 
fully determined by only one time-dependent function. As we shall discover in subsequent sections, there 
exist very interesting dynamical regimes of regular motion depending on more than one frequency and 
characterized by families of quasiperiodic functions, forming s-dimensional tori, with s > 2. 

Such exact regimes of nonlinear mechanical systems, identified by discrete symmetry, are called hushes 
of NNMs. For example, the following two-dimensional bush exists in the FPU-/3 chain with N mod 6 = 0: 



X(t) = {0,2;i(f),X2(t),0,-X2(i),-xi(i)|...|0,xi(t),a;2(i),0,-X2(t),-xi(t)|} 



(20) 



This bush, with multiplicity number m = 6, is determined by two functions, xi{t) and X2{t), satisfying a 
system of two nonlinear autonomous ODEs. 

In contrast to one-dimensional bushes, multidimensional bushes also exist, representing exact 
quasiperiodic dynamical regimes. In the framework of the bush theory, the problem of finding exact non- 
linear regimes in a physical system with discrete symmetry can be solved without any information about 
interparticle interactions. However, the explicit form of the bush dynamical equations essentially depends 
on the interparticle interactions in the considered system. 

Actually, the concept of bushes of NNMs and the group theoretical methods used for their construction 



constitute an extensi on of the theory of complete condensates developed in the papers [Chechin et ami993 



19891 : IChechinl . Il989( | studying phase transitions in solids. Note that the concept of bushes of NNMs does 
not apply only to vibrations of mechanical systems, since one may consider bushes of any other physical 
nature, for example, bushes of spin modes. 

Thus, in sections II, III of this review, we attempt to show that the theory of bushes gives definite and 
quite general answers to the above questions Q1-Q3 in nonlinear systems with discrete symmetries. We 
then proceed to answer question Q4, using the recent analytical and numerical results described in section 
III. 

More specifically, we demonstrate that it is possible to use Poincare-Linstedt series expansions to 
construct exact quasiperiodic solutions with s > 2 frequencies, which belong to the lower part of the linear 
normal mode spectrum, with q = 1,2,3,.... We are thus able to establish the important property that 
the energies Eg excit e d by these low-dimensional, so-called g-tori, are exponentially localized in g-space 
[Christodoulidi et all l2010l |. The relevance of this fact becomes apparent when we study the well-known 
phenomenon of FPU recurrences, which are remarkably persistent in FPU chains (with fixed or periodic 
boundary conditions), at low values of the total energy E. 

We also study the stability of these g-tori as accurately as possible, using methods that do not rely 
on Floquet theory and are quite different than those adopted in section III. In particular, we make use 
of the novel method of the GALI^ indi cators, k = 2, ...,2JV, w hich distinguish rapidly and efficiently 
between chaotic and quasiperiodic orbits Skokos et a/.l . 120071 . l2008l |. As predicted theoretically and verified 
extensively by numerical evidence, in the chaotic case all the GALIs decay exponentially, while in the 
quasiperiodic the first s — 1 GALIs are nearly constant (indicating the dimension s of the associated torus) , 
while all others go to zero by power laws. Thus, by marking the total energy E at which the GALI^, with 
the highest k, start to decrease exponentially, it is possible to determine the destabilization threshold at 
which the g-tori are destroyed, implying the breakdown of FPU recurrences and the eventual equipartition 
of energy among all modes. 



2. Bushes of NNMs in configuration space 



The t heory of bushes of nonl i near n ormal modes was originally developed in jSakhnenko &: Chechinl. 11993 



1994 : IChechin Sz Sakhnenkd . Il998l ]. while its detailed description can be found in Chechin et ali . 120071 ]. 



Low-di mensional bushes in rn echanical sys tems with various kind s of symmetry and st r uctures were stud 
led in [Sakhnenko fc Chechinl . [1993 . 19941 : IChechin fc Sakhnenkd . 1 19981 : IChechin et d\ . l2007l . l2003l . I2OO2I ] . 
In particular, all possible symmetry-determined NNMs, representing o ne-dimensional bushe s , for all A^- 
particle mechanical systems with any of 230 space groups were fou nd in [Chechiri et all 12007 1 . 

The important problem of bush stability was discussed in Chechin et al\ . l2003l . 
first two of these papers being devoted to the vibrational bushes 



2OO2I . l200,^ . the 
in Fermi-Pasta-Ulam chains. Some 
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theorems about bush structu res and the normal forms of their dynamical equations were presented in 
(Chechin Sz Sakhnenkd . Il998l |. Th e general group theoretica l method for the simplification of the bush 



stability analysis was presented in [Chechin Sz Zhukovi . l2006l ] . 



Note that dynamical objects equival ent to the bushes of N NMs have also been discussed for monoatomic 
chains by several other authors [RinkI , l2003l: Poggi Sz Ruffo. 19921- Let us emphasize, however, that the 



group theoretical method s developed in Sakhnenko &: Chechinl . Il993l . 1 19941 : IChechin Sz Sakhnenkd . Il998 



Chechin Sz Zhukov, 2006] can be applied equally w ell not only to monoatomic chains (as illustrated in 



Chechin et al.. . .2002. . .2005i: IChechin fc Zhukovi . l2006l | ) , but also to any other physical system with discrete 



symmetry groups (see Chechin et all 120071 ] ) . 



In this section, we shall demonstrate the power of bush theory methods using three simple mechanical 
models: A nonlinear chain of particles moving in one-dimension under p.b.c, a "square molecule" whose 
particles move in a two-dimensional plane and an octahedral structure of particles moving in three dimen- 
sions. For these nonlinear dynamical systems, we will consider different vibrational regimes and thus deal 
with bushes of vibrational modes. 

Now, in order to properly define bushes of NNMs, one must first determine the parent symmetry group 
of the equations describing the vibrations of a given mechanical system. 

2.1. Parent symmetry group 

A parent symmetry group consists of all transformations which leave the system of dynamical equations 
invariant. In general, such transformations can include space and time variables as well as parameters of the 
system. Let us consider, for example, the dynamical equations ([!]) for the FPU-/? chain with p.b.c. ^ and 
an even number of particles. It is easy to guess some space transformations which do not change Eqs. (j31[6]) 
by considering the equilibrium state of our nonlinear chain. Firstly, this chain is invariant under the action 
of the operator a which shifts it by the lattice spacing a. This operator generates the translational group 

r = {e,a,a^...,a^-^}, a^ = e, (21) 

where e is the identity element and A^ is the order of the cyclic group T. The operator a induces the cyclic 
permutation of all particles of the chain and, therefore, acts on the configuration vector X[t) as follows: 

aX{t) = a{xi(t),X2(t),... ,XN-l{t),XN{t)} = {xN{t),Xi{t),X2{t), . . . ,XM-l{t)}. (22) 

Secondly, the symmetry group of the monoatomic chain contains the inversion i, with respect to the 
center of the chain, which acts on the vector X{t) in the following way: 

lX{t) =l{xi{t),X2{t),... ,XN-l{t),XN{t)} = 

{-XAr(t), -XN-i{t), . . . , -X2{t), -xi{t)}. (23) 

The complete set of symmetry transformations includes also all products a i of the pure translations 
a (/c = 1, 2, . . . , A — 1) with the inversion i and forms the so-called dihedral group D which can be written 
as a direct sum of the two cosets T and T ■ i: 

D = TeT-i. (24) 

This is a non-Abelian group induced by two generators (a and i) with the following generating relations 

a = e, i = e, ia = d~ i. (25) 

When applied to Eqs. (j4]), operators a and i induce the following changes of variables: 

a: Xi{t) ^ X2{t),X2{t) ^ X3{t),...,XN-l{t) ^ XNit),XN{t) ^ xi(t); 
i : Xi{t) ^ -XN{t),X2{t) ^ -XN~l{t),X3{t) ^ -XN-2{t),.-.- 

It is straightforward to check that upon acting on (j4]) with transformation (]26p the system is trans- 
formed to an equivalent form. Moreover, since Eqs. ([4]) are invariant under the actions of a and i, they are 
also invariant with respect to all products of these two operators and, therefore, the dihedral group D is 
indeed a symmetry group of equations (|4]) for a monoatomic chain with arbitrary interparticle interactions. 
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As a consequence, the dihedral group D can be considered as a parent symmetry group for all monoatomic 
nonhnear chains as, for example, the FPU-a chain, whose interparticle interactions are characterized by 
the force 

/(Ax) = Ax + a(Ax)2. (27) 

Of course, since in the case of the FPU-/? chain the force /(Ax) = Ax + /3(Ax)^ is an odd function of 
its argument Ax, the FPU-/3 chain possesses a higher symmetry group than the FPU-a chain. 

Indeed, let us introduce the operator u which changes the signs of all atomic displacements without 
their transposition: 

UX = n{xi(t),X2(t),...,XAr_l(t),XAr(t)} = {-Xi(t), -X2(t), ..., -XAr_i (t), -XAr (t)}. (28) 

It can easily be checked that the operator n generates a transformation of all the variables Xj(t), i = 1..N 
in Eqs. ([I])-©, which leads to an equivalent form of these equations. Therefore, the operator u and all its 
products with elements of the dihedral group D belong to the full symmetry group of the FPU-/3 chain. 

Clearly, the operator u commutes with all the elements of the dihedral group D and we can consider 
the group 

G = D®D-u (29) 

as the parent symmetry group of the FPU-/3 chain. The group G contains twice as m,any elements as the 
dihedral group D and, therefore, possesses a greater number of subgroups. 

On the other hand, every subgroup of the parent group generates a certain bush of NNMs. Therefore, 
there exists a greater number of bushes for the FPU-/3 chain (and for any other chain with odd force 
of interparticle interactions) compared with those of the FPU-a chain (and for all arbitrary nonlinear 
monoatomic chains, as well). 

Finally, let us note that it is sufficient for our purposes to define any symmetry group by listing only 
its generators which we denote by square brackets, for example, we write T[a], D[a,i], G[a,i,ii]. 

2.2. Subgroups of the parent group and bushes of NNMs 

Let us consider now a specific configuration vector X^^'{t) [see Eq. ([7])], which determines a displacement 
pattern at time t, and let us act on it successively by the operators g that correspond to all the elements 
of a parent group Gq. The full set Gj of elements of the group Go, under which X^^>{t) turns out to be 
invariant, generates a certain subgroup of the group Gq (Gj C Gq) [Note that this subgroup can be trivial, 
i.e. it can consist of only one symmetry element e. In this case, a chosen configuration vector X^^'(t) turns 
out to be of general form, since there are no connections between the displacements of different particles of 
the chain]. We then call X^^'{t) invariant under the action of the subgroup Gj of the parent group Gq and 
use it to determine the bush of nonlinear normal modes corresponding to the subgroup Gj of the group 
Gq. 

Therefore, in the framework of the above approach, one must find all the subgroups of the parent group 
Go to obtain all the bushes of NN Ms of different types,. This can be done by standard group theoretical 



methods. In [Chechin et all 120051 ] a simple crystallographic technique was developed for singling out all 
the subg roups of the pa.rent group of any rnqnoato mic chain, following the approach of a more general 
method Chechin et all Il993l . Il989l : IChechinl . 119891 ] . We shall not discuss here how to find subgroups of 



the parent groups; rather, we will demonstrate how one can obtain bushes of NNMs if the subgroups are 
already known. 

Let us consider the subgroups Gj of the dihedral group D. Each group Gj contains its own translational 
subgroup Tj C T, where T is the above discussed full translational group pij) . If N is divisible by 4 (for 
example, we consider below the case A^ = 12) there exists a subgroup T4 = [a^] of the group T = [a]. Note 
that, in square brackets, we write down only the generators of the considered group, while the complete 
set of group elements is written in curly brackets, as, for example, in Eq. (I2ip . 

If a vibrational state of the chain possesses the symmetry group T4 = [a^\ = {e, a^, a^, . . . , a^~^}, the 
displacements of the atoms, which are at a distance 4a from each other in the equilibrium state, turn out 
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to be equal, since the operator a^ leaves the vector X{t) mvariant. For example, for the case A^ = 12, the 
operator a'^ permutes the coordinates of X = (xi, X2, ..., 2;i2) taken in quadruplets (xj, Xj+i, Xj+2, Xj+a), i = 
1,5,9, while from equation d^X(t) = X{t) one deduces Xi = Xi+4,i = 1,2,3,4. Thus, the vector X{t) 
contains 3 times the quadruplets xi,X2,X3,X4, where Xj(t) (i = 1,2,3,4) are arbitrary functions of time 
and can be written as follows 

X{t) = { Xi{t),X2{t),X3{t),X4^{t) I Xi{t),X2{t),X3{t),X4{t) \ Xi{t), X2{t), X3{t), X^it) }. (30) 

In other words, the complete set of atomic displacements can be divided into N/4: (in our case, N/A = 3) 
identical subsets, which are called "extended primitive cells" (EPC). In the bush (pOl) . the EPC contains 
four atoms, and the vibrational state of the whole chain is described by three such EPC. Thus, the EPC 
for the vibrational state with the symmetry group T4 = [a'^] has size equal to 4a, which is four times larger 
than the primitive cell of the chain in the equilibrium state. 

It is essential that some symmetry elements of the dihedral group D disappear as a result of the 
symmetry reduction D = [a,i] — )• T4 = [a^]. There are four other subgroups of the dihedral group D, 
corresponding to the same translational subgroup T4 = [a^]: 

[a ,i], [a ,ai], [a ,a i], [a ,a i]. (31) 

Each of these subgroups possesses two generators, namely d^ and an inversion element d i {k = 0, 1, 2, 3). 
Note that the d i differ from each other by the position of the center of inversion. 

Subgroups [d'^,d'^i] with A; > 3 are equivalent to those listed in (f3T]l . since the second generator d'^i 
can be multiplied from the left by d~^, representing the inverse element with respect to the generator d^. 
Thus, there exist only five subgroups of the dihedral group (with N mod 4 = 0) constructed on the basis 
of the translational group T4 = [d^] , namely, this group and the four groups from the list (|3T]l . 

Now, let us examine the bushes corresponding to the subgroups (f3T]l . The subgroup [d^,i] consists of 
the following six elements: 

^ ^4 ^g ^ ^4^ „8„ _ ^^4 /„„\ 

e, a , a , 1, a 1, a 1 = la . [62) 

The invariance of X(t) with respect to this group can be written as follows: 

a^X{t) = X{t), iX{t) = X{t), (33) 

while the invariance of the vector X{t) under the action of the group generators [d^ and i] guarantees its 
invariance under all elements of this group. 

As explained above, the equation a^X{t) = X{t) is satisfied by the vector X{t), see (pO]l . while iX{t) = 
X{t) also holds, from which we obtain the following relations for each of the three EPS: xi{t) = —Xi{t), 

X2{t) = -X3{t). 

Therefore, for N = 12, the invariant vector X{t) of the group [d'*, i] can be written in the form 

X{t) = {xi{t),X2{t),-X2{t),-Xi{t)\xi{t),X2it),-X2{t),-Xi{t)\xi{t),X2{t),-X2{t),-Xi{t)}, (34) 

where xi{t) and X2{t) are arbitrary functions of time. 

Thus, the subgroup [a^,i] of the dihedral group D generates a two-dimensional bush of NNMs. The 
explicit form of the differential equations governing the two variables xi{t) and X2{t) can now be obtained 
by substitution of the ansatz ^^ into the FPU-/3 dynamical equations (jl][6]). We shall, hereafter, denote 
the bush (j34p in the form 

B[a'^,i] = \xi,X2,-X2,-xi\, (35) 

showing the atomic displacements in only one EPS and omitting the argument t of variables xi{t), X2(t). 
Proceeding in a similar manner, we obtain bushes of NNMs for three other groups from the list (I31J1 : 

B[a ,ai] = |0, x,0,— x|, (36) 

i?[d^,d^i] = |xi,-xi,X2,-X2|, (37) 
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B[a 



4 a^] 



\x, 0, —X, 



(38) 



Let us comment on these results: The bushes B[a'^,di] and B[a'^,a'^i] turn out to be one-dimensional 
in the sense that they represent simple periodic solutions corresponding to NNMs. These modes are, in 
fact, equivalent since they differ from each other only by the numbering of the particles on the chain. In 
group theoretical terms, this equivalence is a consequence of the fact that Gi = [a^,ai] and G2 = [a'^,a^i] 
prove to be conjugate |j subgroups in the parent group D = [a,i]. Both two-dimensional bushes ([35]) and 
()37p are also equivalent to each other for the same reason. 

All displacement patterns (I35II38P are particular cases of the general pattern of the bush Bla'^] = 
\xi,X2,X3,xA. One may, therefore, ask what will happen during time evolution of other particular cases 
of i?[a^]. For example, if one chooses for the FPU /3-chain with A^ mod 4 = an initial one-dimensional 
pattern of the form X{0) = {3xo, —xq, —xo\3xo, —xq, —xo\---} or, for N mod 5 = 0, an initial state of 
the type X{0) = {—xo,—xo,4xq,—xo,—xo\—xq,—xo,4:Xo,—xq,—xo\...} one easily finds that there are no 
subgroups of the parent group D which produce these displacement patterns. As a consequence, when 
solving the equations of motion of the FPU-/3 chain in these cases to obtain X{t), the above structures will 
not be conserved in time. 

More generally, we conclude that for sufficiently large EPS there are not enough symmetry elements to 
give rise to NNMs, since the bushes of the corresponding displacement patterns are multidimensional. For 
this reason, there exists only a very specific number of bushes for any fixed dimension beyond the NNMs! 



2.3. Bushes of NNMs for arbitrary interparticle interactions 

As was proved in [Chechin et all |2005| | for any nonlinear monoatomic chain with arbitrary interparticle 
interactions, there exist only three symmetry-determined NNMs: 



B[a 



a2 



-X 



, (vr-mode), B[a^ 



x,0, 



-X 



B[a ,ai 



0,2;, 0, 



(39) 



(each NNM corresponds to a certain subgroup of the dihedral group D). These modes can easily be excited 
in the FPU-a and FPU-/3 chains, if we solve numerically the corresponding dynamical equations with 
appropriate initial conditions. 

On the other hand, it is already known that the FPU-/3 chain possesses a higher symmetry group G, 
since changing signs of all the displacements by the operator u [see Eq. (j28p ] leaves the FPU-/3 dynamical 
equations invariant. Therefore, all subgroups of the group D also turn out to be subgroups of the group 
G, while G possesses some additional subgroups whose elements can contain the operator -u. 

As a consequence, there exist more bushes of NNMs for the FPU-/? chain (and for any other chain 
with odd force of interparticle interactions) than for the FPU-a chain (or any other chain whose inter- 

{ )article interactions is not described by odd function /(Ax)). These additional bushes were presented in 
Chechin et al\ . |2005| | in the following form: 



B[a^,iu] 



X, — 2x, X 



B\a?uA 






, B\pru^ai 



\j • tXj • tij s \j m. tXj • ■/> 



(40) 



and were also obtained in Rinkl . l2003l | by a different method 



2.4. Bushes of NNMs for a square molecule 

In our model, a square molecule is represented by a mechanical system whose equilibrium state is depicted 
in Fig. [21 Four atoms of this molecule are shown as filled circles at the vertices of the square. The number 
of every atom and its (x, y) coordinates are as given in the figure. 

Let us suppose that atoms can oscillate about their equilibrium positions only in the (x, y) plane and, 
therefore, eight degrees of freedom are needed to describe this mechanical system. Furthermore, we will 



Two subgroups Gi and G2 of the same group G are called conjugate to each other, if there exists at least one element qq £ G 
which converts Gi into G2 via the transformation G2 = Qq Giqq. [note that in the case of bushes ([36)l and (|38|l . go = a]. 
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Fig. 2. Mechanical model of a square molecule. 



not consider any specific type of interparticle interactions and treat bushes of NNMs as purely geometrical 
objects. The equilibrium configurations of our molecule possess the symmetry group denoted by C^^ in the 
Schoenflies notation, while, at equilibrium, the molecule depicted in Fig. [2] is invariant under the action of 
the following transformations: 

Rotations through the angles 0°, 90°, 180°, 270° about the z axis orthogonal to the plane of Fig. [2] and 
passing through the center of the square. We will denote these rotations by gi, §2, 53, Qa, respectively. 
Reflections in four mirror planes orthogonal to the plane of Fig. [2] and passing through the z axis. (We 
depict these mirror planes by bold lines). Two of them are "coordinate" planes ((75, gj) and the other two 
are "diagonal" planes {g^, g^). 



Analytically, the above mentioned symmetry elements can be defined as follows: 



51(2;, y) = {x, y), 52(2;, y) = {-y, x) 
95{x, y) = (-X, y), gQ{x, y) = {-y, 



93{x, y) = (-X, -y), 54(2;, y) = {y, -x) 
-x), g7{x,y) = {x,-y), gs{x , y) = {y , x) 



(41) 



Thus, the symmetry group of the square molecule (C^v) contains 8 elements gi, ■ ■ ■ ,g8, determined 
by Eq. (J4T]) . This group is non-Abelian since, for example, 52 • ffs = 95, while g% ■ g2 = gj- According to 
Lagrange's theorem, the order of any subgroup is a divisor of the order of the full group. Therefore, for the 
case of the group G = C^y with the order m = \\G\\ = 8, there exist only subgroups Gj with order equal 
to m = 1, 2, 4, 8 given by the following list: 



m 
m 



m 



m 



1 : Gi = {gi} = Ci 
2: G2 = {51,53} = C2 

^3 = {51,55} and G3 

Gi = {51,56} and G'4 
4 : G5 ={51,52,53,54} = 

Gq = {51,53,55,57} = 

Gj = {51,53,56,58} = ^2v 
8: Gg = {51,52,53,54,55,56,57,58} 



{51,57} 
{51,58} 






r-ic 



(42) 



G 



4« 



Note that throughout this paper, we use Schoenflies notation for the point symmetry groups and 
distinguish between coordinate and diagonal settings of some subgroups using the superscripts 
"d", respectively. 
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Let us suppose that the equihbrium state of our molecule is stable under arbitrary infinitesimal dis- 
placements of the atoms in the xy-plane. Moreover, we will also assume that this state is isolated in the 
sense that within a finite radius neighborhood around it there are no other equilibrium states. 

Next, we consider the planar vibrations, i.e. vibrations of the molecule in the plane of its equilibrium 
configuration. Let us excite a vibrational regime of our molecule by displacing the atoms from their equi- 
librium positions in a specific manner. As a result of such displacements, the initial configuration of the 
molecule will have a well defined symmetry which is described by one of the subgroups of the group G4t, 
listed in (02]). 

Indeed, the first configuration in Fig. [3] with atoms displaced arbitrarily corresponds to the symmetry 
group Gi = Ci. In this figure, we depict by arrows the atomic displacements and by thin lines the resulting 
instantaneous configurations of the molecule. 

For our analysis, it is essential that the symmetry of the instantaneous configuration of the mechanical 
system be preserved during the vibrational motions. More precisely, if an atomic pattern at time to has a 
symmetry element g, this element cannot disappear spontaneously at any t > tQ. This proposition, which 
we call the "symmetry preservation theorem", can be proved rigorously by examining the Hamiltonian 
equations of motion. 

Let us note, however, that spontaneous lowering of the symmetry can occur when the considered 
dynamical regime loses its stability leading most frequently to the appearance of another bush of higher 
dimension. This important phenomenon, which may be regarded as the dynamical analogue of a phase 
transition, will be discussed in Section III. 

Of course, at the times when all atoms pass through their equilibrium posit i ons, t he molecule config- 



uration is again a square. Nevertheless, it can be proved Sakhnenko &: Chechinl . Il993l | that these isolated 
instances do not affect such general dynamical properties as the selection rules for excitation transfer 
between modes of different symmetry. 

Thus, all possible vibrational regimes of the square molecule can be classified according to 8 subgroups 
of the group G = C^y. The different configurations of the vibrating molecule, therefore, as depicted in 
Fig. [3l are the following: a pulsating square (Gg = C^v), a rotating and pulsating square (G5 = C4), a 
rectangle {Gq = G^^), a rhombus (G7 = G2^), a trapezoid (G3 = G^), a deltoid (G4 = Gf ), a parallelogram 
(G2 = G2) or an arbitrary quadrangle (Gi = Gi). All these configurations vary in size as time progresses, 
but the type of the corresponding quadrangle does not change. 

Note that the different types of the above vibrational regimes of the molecule are described by different 
numbers of degrees of freedom. Let us discuss this important fact more carefully: 

Each of the eight types of vibrational regimes in Fig. [3] corresponds to a certain bush of NNMs. For 
example, the dynamical regime representing a pulsating square (G = Gi^) can be characterized by only one 
degree of freedom. The edge of the square or the displacement of a certain atom from its equilibrium position 
along the corresponding diagonal can be used to represent this degree of freedom. Such a vibrational regime 
is described by a one- dimensional bush consisting only of the so-called "breathing" nonlinear normal mode. 

The rhombus- type vibration (G7 = G|„) and rectangle- type vibration (Gg = G^^,) are characterized by 
two degrees of freedom. The length of the diagonals in the former case and the length of the adjacent edges 
in the latter can be chosen as these degrees of freedom. Both of these vibrational regimes are described 
by two-dimensional bushes of modes. Similarly, one can check that a trapezoid-type vibration (G = G^) 
corresponds to a three-dimensional bush, the deltoid-type vibration to a four- dimensional bush and the 
vibration with arbitrary quadrangle is represented by a five- dimensional bush of vibrational modes. Note 
that vibrations with G5 = G4 group are described by a two-dimensional bush consisting of one rotating 
and one pulsating mode. 

Let us also emphasize that the above-discussed classification of vibrational regimes does not depend on 
the interparticle interactions characterizing the mechanical system. On the other hand, the time-evolution 
of the degrees of freedom corresponding to a certain bush does depend on these interactions in an essential 
way, as we discuss in more detail below. 
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1) Ci 




5) Q 




2) Co 




6) a 
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2v 




3) C-, 




7) ci 




4) C', 




a 



4v 




Fig. 3. DifTerent vibrational regimes (bushes of NNMs) of the square molecule. 
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2.5. Bushes of NNMs for a simple octahedral structure 

Finally, let us describe the occurrence of bushes in a S-dimensional mechanical system consisting of six 
mass points (called particles or atoms) whose interactions are described by a pair isotropic potential u(r), 
r being the distance between two particles. We suppose that, at equilibrium, these particles form a regular 
octahedron with edge qq as depicted in Fig. HI Let us introduce a Cartesian coordinate system. Four 
particles of the above octahedron lie in the x, y plane and form a square with edge uq. Two other particles 
lie on the z-axis and are called "top particle" and "bottom particle" with respect to the direction of the 
z-axis. The distance between these particles is also equal to gq. 




Fig. 4. Octahedral mechanical system. 



At equilibrium, this syste m possesses the poi nt symmetry group Oh and all its possible bushes of 
vibrational modes are listed in Chechin et aLl . l2003l |. We may suppose that u{r) is the well-known Lennard- 
Jones potential: 

^(^) = ^ - ^' (43) 

as we have done in Chechin et ali |2003|, where we studied the stability of the bushes in the octahedral 
molecule. In the present section, however, we shall consider u{r) as an arbitrary pair isotropic potential. 
The potential energy of our system in its vibrational state is written in the form 



UiX) = Y, <nj 



uv' 



(44) 



1,3 



where rij is the distance between the ith and jth particles. The A^-dimensional vector X = 
{xi(t),X2(t), . . . ,XAr(i)} in Eq. ^^ determines the displacements of all particles at an arbitrary chosen 
instant t. Here A^ = 18 is the total number of degrees of freedom of the octahedral mechanical system. 

The configuration vector X{t) = {xi{t) , X2{t) , . . . ,xis{t)} determines all the dynamical variables Xi{t) 
as follows: The first three components of this vector correspond to the x-, y- and z-displacements of particle 
1, the next three components correspond to x-, y-, z-displacements of particle 2, etc., as numbered in Fig. 
m The dynamics of this system is described by Newton's equations 

dU 

Xi = --—, ^ = l,2, ...,18 

OXi 



with the masses of all particles equal to unity. 
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Each bush of modes corresponds to a certain subgroup Gj of the symmetry group Go = Oh of the 
mechanical system at equiUbrium. This means that the vibrational state described by the bush possesses 
the symmetry group Gj C Go and, therefore, there exist certain restrictions on the dynamical variables 
Xi{t). As a result, the number m of degrees of freedom corresponding to the given bush (i.e. its dimension), 
is smaller than the total number of dynamical variables of the system. 

Indeed, at any time t, the atom configuration in a vibrational state represents a polyhedron char- 
acterized by the symmetry group Gj C Oh- Thus, instead of the old variables Xi(t), we introduce m 
new variables yj{t) (j = l,...,m), which compl etely determine th i s poly hedron. Next, we employ the 
well-known Lagrange method (see, for example. Landau Sz Lifshita . Il976l |) for obtaining the dynamical 



equations in terms of the new variables yj{t). 

For this purpose, we introduce the Lagrange function L = T — U, where T is the kinetic and U is 
the potential energy expressed as functions of new variables yj and new momenta yj (j = 1 , . . . , m) . Then 
Euler-Lagrange equations in the new variables read 

where m is the dimension of the considered bush. 

We will not present here the tedious derivation of such equations for the octahedral molecule (see 
e.g. Chechin et aU 120031 ] . where dynamical equations for bushes are obtained by means of an appropriate 



MAPLE-program) . Rather, we will give the final results for the bushes B[0/i], B [1)4/1] and B[C4^]. Here, in 
the square brackets next to the bush symbol B, we give the symmetry group of the corresponding bush 
of vibrational modes. The symmetry groups of the above bushes satisfy the following group-subgroup 
relations: 

C4. C D^h C Oh- (46) 

The bushes B[0/i], B[L)4/j], and B[G4^] have one, two and three dimensions, respectively, while their 
geometrical features can be immediately seen from the symmetry groups that correspond to them. Indeed, 
the one-dimensional bush B[0/i] consists of only one ("breathing") mode: The appropriate nonlinear dy- 
namical regime describes the evolution of a regular octahedron whose edge a = a{t) periodically changes 
in time. 

The two-dimensional bush B [1)4/1] describes a dynamical regime with two degrees of freedom. The 
symmetry group G = D^h of this bush contains the four- fold axis symmetry coinciding with the z coordinate 
axis and the mirror plane coinciding with the x, y plane. This symmetry group restricts essentially the shape 
of the polyhedron describing our mechanical system in the vibrational state. Indeed, the presence of the 
four-fold axis demands that the quadrangle in the x,y plane be a square. For the same reason, the four 
edges connecting the particles in the x, y plane (vertices of the above square) with the top particle lying 
on the z-axis must have the same length, which we denote by b(t). 

Similarly, let the length of the edges connecting the bottom particle on the z-axis with any of the four 
particles in the x, y plane be denoted by c{t). In the case of the bush B[L)4/i], b{t) = c{t) for any t, because of 
the presence of the horizontal mirror plane in the group G = D^h- However, for the three-dimensional bush 
B[C4„], this mirror plane is absent and, therefore, b{t) ^ c{t). We illustrate the instantaneous configuration 
of our mechanical system vibrating according to the bush B[G4i,] in Fig. [5j 

Let us also introduce the heights, hi{t) and h2{t) of the perpendicular distances, respectively, from the 
top and bottom vertices of our polyhedron onto the x, y plane. We can now write the dynamical equations 
of the above bushes in terms of the purely geometrical variables a{t), b{t), c{t), hi{t) and /i2(i). 

Choosing a{t) and h{t) = hi{t) = /i2(t) as suitable variables for describing the two-dimensional bush 
B[Z)4/i] and a{t),hi{t) and /i2(t) as appropriate for describing the three-dimensional bush B[G4^], we can 
write the potential energy for our bushes of vibrational modes as follows: 

B[Oh] : Uia) = 12u(a) + 3ii(\/2a), 

B[D4h] : U{a,h) = Au{a) + 2u{V2a) + 8u ( ^Jh^ + ^ J +n(2/i), 

B[C4^] : U{a, /ii, /12) = 4u(a) + 2u{V2a) + Au{b) + Au{c) + u{hi + /i2). 
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Fig. 5. The distorted octahedron illustrating positions of particles at fixed instant for the bush B[C4t,] 



where b 



4 + ijhi-^h2) , c 



— 4- 



{lh2 



-M 



Then, with the aid of Lagrange's equations, we obtain the fohowing dynamical equations for the above 
bushes of vibrational modes: 



B[Oh] : 

a = 

B[D^h] 
a = 
h = 

B[C4v] '■ 

a = 
hi- 
h2-- 



-4u'{a) - V2u'{V2a); 

-2u'{a) - ^/2u'{^/2a) - 2n'(6)f , 
-4n'(6)f -u'(2/i); 

-2u'{a) - V2u'{V2a) - u'{b)^ - ^'(c)! , 
u'{b)^-hi^-u'ih, + h2), 
u'{hi + h2). 



(47) 



.U'(c)^^i^^^ 



Thus, we derive dynamical equations for our bushes of vibrational modes in terms of variables having 
an explicit geometrical nature. Each bush describes a certain nonlinear dynamical regime corresponding to 
a vibrational state of the considered mechanical system, so that at any fixed time the configuration of this 
system is represented by a definite polyhedron characterized by the symmetry group G of the given bush. 

The above dynamical equations for the bushes B[0/i], B[L'4/i] and B[C4„] are exact, but rather compli- 
cated to solve. One may consider instead some approximate equations which can be obtained by expanding 
the potential energy in Taylor series (near the equilibrium state) keeping only the lowest order terms in 
the expansion. 

Thus, using the first few leading terms in such decompositions, it turns out that bushes belonging 
to systems of different physical nature, with different symmetry groups and different structures can pro- 
duce equivalent dynamical equations. More precisely, the equations for many bushes can be written in 
the same form, and this fact leads to the idea of the "classes of dynamical universality" mentioned in 
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Returning to the dynamical equations (j47l) for the bushes B[0/i], B [1)4/1], B[C4t,], we see that every 
bush can be considered as a reduced dynamical system whose dimension [m) is less than the full dimension 
of the considered mechanical system. If the original system is Hamiltonian, the reduced system can also 
be proved to be Hamiltonian. 

3. Bushes of NNMs in modal space and stability analysis 

Vectors X{t), corresponding to bushes of NNMs, are determined in the configuration space Jft . If we 
introduce in this space a basis set of vectors {(pi,(p2-, ■■■■,'^n}, the dynamical regime of our mechanical 
system can be represented by a linear combination of the vectors (pj with time dependent coefficients i'j{t) 
as: 

N 

X{t) = Y,^,{t)0, (48) 

i=i 

where the functions ^^{t) entering this decomposition represent new dynamical variables. 

Every term in the sum ()48p has the form of a NNM, whose basis vector (pj determines a displacement 
pattern, while the functions yj{t) determine the time evolution of the atomic displacements. Because of 
this interpretation, we can consider a given dynamical regime X[t) as a hush of NNMs. In fact, as will be 
shown later one often speaks about root modes and secondary modes of a given bush. 

Note that each NNM i'j{t)(pj from ()48p is not, in general, a solution of the dynamical equations of the 
considered mechanical system, while a specific linear combination of a number of these modes can represent 
such a solution and, thus describe an exact dynamical regime. Sometimes, for brevity, we will use the word 
"mode" not only for the whole term Vj{t)(pj, but also for the time dependent function ^^{t). 

3.1. Normal coordinates and normal modes 

Normal coordinates and normal modes are commonly introduced in the theory of small vibrations of 
multiparticle systems (see, for example. Landau &: Lifshita . Il976l |) in the framework of the harmonic 



approximation, where they represent exact solutions of the corresponding dynamical equations. 

In this approximation, the potential energy U{X) is decomposed into a Taylor series with respect 
to the atomic displacements Xj(t) from their equilibrium positions and all terms whose orders are higher 
than 2 are neglected (because the displacements are supposed to be sufficiently small). As a result, the 
dynamical equations 

dU 

miXi = -j—, {i = l,...,N) (49) 

OXi 

turn out to be linear differential equations, with constant coefficients. 
Each normal mode is a particular solution to Eq. (I49p of the form 

X{t) = ccos{ujt + (fo), (50) 

where the A^-dimensional constant vector c = {ci,C2, . . . , cat} and the constant phase (pQ determine the 
initial displacements of all particles from their equilibrium state and cj is the frequency of vibration. 

Since U{X) has a quadratic form, substituting ([50]) into Eq. (i49|) and dividing the resulting equations 
by cos{ujt + fo) , we reduce the problem of finding the normal modes to the algebraic problem of evaluating 
the eigenvalues and eigenvectors of the matrix /C = ||A;jj|| with coefficients 

i,j = l..N. (51) 

x=o 

Each eigenvalue of this matrix is the square of the frequency, w?, of the jth normal mode, while its 
eigenvector Cj represents the pattern of this normal mode and is called the normal coordinate. 



^'' dxjdxi 
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Being real and symmetric, the N x N matrix /C possesses A^ eigenvectors Cj {j = 1, . . . ,N) and N 
eigenvalues a;|. The complete collection of these vectors, i.e. the N normal coordinates, can be used as the 
basis of the configuration space, hence we may write 

N 

Xit) = ^f,,{t)c„ (52) 

«=i 

where X{t) = {xi(t),X2{t), . . . , xi\i{t)}, while /Xj(t) are new dynamical variables which we introduce instead 
of the old variables Xi{t) (i = 1, . . . , A^). 

If the transformation to normal coordinates is used in the absence of degeneracies, the corresponding 
system of linear ODEs leads to a set of uncoupled harmonic oscillators 

il,{t)+oj]fij{t) = 0, j = l,...,N 

with the well-known solution 

IJ,j{t) = aj cos {cojt + ifoj) , (53) 

where Uj and ipoj are arbitrary constants. 

Note the distinction we make between a normal coordinate, represented by the eigenvector Cj, and a nor- 
mal mode, referring to the product of the vector Cj and the time-periodic function /Xj(t) = cos {cojt + (foj)- 

3.2. Symmetry -adapted coordinates as basis of the configuration space 

In general, the normal coordinates (normal modes) depend on the interparticle interactions in the mechan- 
ical system. Indeed, the concrete form of the potential energy is needed for their construction according to 
the above-described procedure. Because of this fact, it is more convenient to use symmetry -adapted coor- 
dinates, instead of normal coordinates. The former are the basis vectors of the irreducible representations 
of a certain parent group. More importantly, the symmetry-adapted coordinates do not depend on the 
interparticle interactions and can be obtained by the appropriate group theoretical methods. 

Note that, for simple mechanical systems, the normal and the symmetry-adapted coordinates may 
coincide with each other. This coincidence is a consequence of the fact that each irreducible representation 
of the parent group is contained only once in the decomposition of the reducible vibrational representation 
into its irreducible parts. 

We complete this part of our study by recalling in Appendix A some notions from the theory of matrix 
representations of symmetry groups. We also discuss in Appendix A the theory of irreducible representations 
and describe how they can prove especially helpful in analyzing important problems of crystal vibrations 
in Solid State Physics. 

This concludes our brief outline of the group theoretical concepts needed to find bushes of 
vi brational modes for arbitra r y mechanical systems. The interested reader will find more details 
in 



Sakhnenko fc Chechinl . Il99,i IChechin fc Sakhnenkol . Il998l : IChechin et 'ail . 120001 : IChechinl . Il989l | . Now,it 



is time to discover how this approach of discrete symmetries can help us analyze the stability properties 
of the periodic and quasiperiodic solutions belonging to these bushes. 

3.3. Linear stability analysis of bushes of NNMs 

Let us begin by considering individual NNMs, representing one-dimensional bushes of the FPU-/3 Hamil- 
tonian. 

1. Example 1: Stability of the 7r-mode in the FPU-/? chain 

The question of local stability near NNMs of the form Xj may be answered by studying the linearized 
equations about them, setting Xj = Xj + yj and keeping up to linear terms in yj. For example, in the case 
of the NNM ([8]) these equations become (/3 = 1): 

jjj = (1 + A8x^){yj-i - 2yj + y^+i), j = 1, . . . , AT (54) 
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where x is the solution of ([9|), see also equations pO|) - (fT3|) . and yi = un+i-, Ho = UN, due to p.b.c. 

According to Floquet theory, linear stability analysis of periodic solutions is performed by studying the 
eigenvalues of the monodromy matrix M(T), whose N columns are the fundamental solutions of Eq. (j54p 
evaluated at t = T (the period of the NNM), with M(0) equal to the N x N identity matrix. However, in 
the case of such simple modes as ([8]) , it is much easier to diagonalize first the variational equations into A^ 
uncoupled Lame equations 

z,{t) + 4{l + A8x^)sin^(^yj{t) = 0, j = l,...,N (55) 

where the Zj variations are simple linear combinations of the y^'s. Changing variables to u = At, where 
A^ = 4/(1 — 2k'^) and using ()10p -()13 p this equation takes the form 

Zj{u) + {l + 22k'^ -24.k'^sn'^{u,k'^))sm^{^)zj{u) =0, j = l,...,N (56) 

where primes denote differentiation with respect to u. Equation (j56p is an example of Hill's equation 

z"{u) + Q{u)z{u) = (57) 

where Q{u) is a T-periodic function {Q{u) = Q{u + T)). According to Floquet theory, its solutions are 
bounded or unbounded depending on whether the eigenvalues of its monodromy matrix lie on the unit 
circle or not. The first variation Zj{u) to become unbounded as k > kc determines the total energy of the 
lattice 

at which this particular NNM becomes unstable for < /c < -72 

2. Example 2: Stability of the 7r-mode in the FPU-a chain 

Let us consider the stability of the one-dimensional bush B[d?,i\, representing the vr-mode, in the 
FPU-a chain with N = 4 particles under p.b.c. 
The equations for this nonlinear chain read 

Xl = f{x2 - Xl) - f{xi - X4), 

X2 = f{x3-X2)- f{x2-Xi), .^g. 

X3 = f{Xi-X3) - f{x3-X2), 

XA = f{xi - X4) - /(X4 - X3), 

where f{x) = x + ax"^ . Clearly, the vr-mode represents here the exact periodic vibrational regime X = 
{A{t) , — A{t) , A{t) , — A{t)} , with A{t) = Cocos(2t). According to the conventional prescription, we should 
linearize the dynamical system (|59p in the neighborhood of the given bush and study the resulting system 
of ODEs. To achieve this, we write 

X{t) = C{t) + 5{t), (60) 

where (7 = {A{t) , - A{t) , A{t) , - A{t)] is our NNM, while ^(t) = {5i{t),52{t),53{t),5i{t)] is an infinitesimal 
vector. Substituting ([60]) into Eqs. ([59]) with /(x) = x + x^ (i.e. a = 1) and neglecting all terms nonlinear 
in 5j(t), we obtain for this FPU-a model the system of linearized equations: 

5 = J{t)l (61) 

where J{t) is the Jacobi matrix of (I59p for the FPU-a chain calculated by the substitution of the vector 
X = C = {A(t), —A{t),A(t), —A{t)}. This matrix can be expressed as follows: 

J{t) = C + 4A{t)M, (62) 
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(63) 



are time-independent symmetric matrices. 

It easy to verify that C and Ai commute with each other; therefore, there exists a time-independent 
orthogonal matrix S transforming both £ and Ai into diagonal form: SCS = Cdia-, SM.S = M.dia (where 
S is the transpose of S). In turn, this means that the Jacobi matrix J{t) can be diagonalized at any time 
t by one and the same time-independent matrix S. Therefore, our linearized system (j6ip for the considered 
bush B[a^,i] can be decomposed into four independent differential equations. 

Observe that, due to similar reasons, we were able to reduce ()54p into a system of N uncoupled ODEs, 
see (j55p . for the vr— mode of the FPU /3— chain in Example 1. 

The above matrix S can be obtained with the aid of the theory of irreducible representations of the 
symmetry group G, for the case G = D2 considered here. The general method for obtaining the matrix S 
uses the basis vectors of irreducible representations of the group G, constructed in the mechanical space 
of the considered dynamical system. In the simple example of a monoatomic chain with A^ = 4 particles, 
this method leads to the following result 





/I 1 1 1 
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1 1-1-1 


2 


1-1 1-1 




I 1-1-1 1 



(64) 



The rows of the matrix S in (j64p are simply the characters of four one-dimensional irreducible representa- 
tions (irreps) - Fi, r2, Ta, r4 - of the Abelian group D2, because each of these irreps is contained once in 
the decomposition of the natural representation of the group G = D2 (see also Appendix A). Introducing 
new variables Y = {yi,y2;?/3;?/4} by the transformation Y = S6, with S as in (j64p . we arrive at the 
com,plete splitting of the linearized equations ()6ip for the FPU-a model: 

(65a) 
(65b) 
(65c) 
(65d) 



?/i = 0, 

m = -2[1 + 4A{t)]y2, 

m = -4y3, 

i/4 = -2[1 - 4A(i)]2/4, 



where A{t) = Gq cos(2t). 

With the aid of Eqs. (I65p one can now determine precisely the threshold value for Co, leading to the 
loss of stability of the one-dimensional bush B[a^,i]. Indeed, according to Eqs. (f65]l . the variables yj{t) 
{j = 1,2,3,4) are independent from each other, and we can solve for them separately: Eq. (I65ap for yi{t) 
describes the uniform motion of the center of mass of our chain, since it follows from the equations Y = S6 
that yi{t) = {Si{t) + 62{t) + 63{t)-\-64{t))/2. Therefore, considering vibrational regimes only, we may assume 
yi{t) = 0. 

If, in the solution of the system (i65]l . yj{t) = for j = 1,2 but not j = 3, 5 = SY would imply 
(since S is orthogonal and S~^ = S): 5{t) = {y3(t), — y3(t),y3(t), — y3(t)}, where y3{t) ~ cos(2t). This leads 
to variations "along" the bush X{t) = Co{cos(2t), — cos(2t),cos(2t), — cos(2t)}, which do not affect the 
stability of the motion. 

On the other hand, if y2{t) and y4(t) are nonzero, with A{t) = Gocos(2t), Eq. (I65bp and (I65dp imply 
ijj + [2 + 8Cn cos(2t)]ii/.,- = 0,7 = 2,4 and can thus be transformed to the standard form of Mathieu's 
equation Abramowitz Sz StegurJ . Il965l ]. Therefore, the stability threshold of the considered bush B[a^,i] 
for A^ = 4 can be determined directly from the well-known diagrams of stable and unstable motion of the 
Mathieu equation. In this way we find that the lowest critical value Cc > of the amplitude Gq at which 
the given bush loses its stability is Gc ~ 0.303. 
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Let us now consider the FPU-a chain with A^ = 6 particles. Linearizing the dynamical equations of 
the FPU-a chain with A^ = 6 in the vicinity of the bush B[d^,i] (vr-mode), we also obtain a Jacobi matrix 
of the form (j62p . However, unlike the case A^ = 4, matrices C and Ai, do not commute with each other. 
As a consequence, we cannot diagonalize both matrices C and Ai, i.e. with one and the same orthogonal 

matrix S. Thus, it is impossible to diagonalize the Jacobi matrix J^{t) in equation 5 = Cf{t)5 for all time t. 
This means that there is no matrix S that completely decouples the linearized system for the bush B [o^ , i] 
for the chain with A'^ = 6 particles. 

This difference between the cases A^ = 4 and A^ = 6 (generally, for A'^ > 4) can be explained 
group-theoretically as follows: The group G = [a^,i] of the considered bush leads to different groups 
for the cases A^ = 4 and A^ = 6. Indeed, for A^ = 4 [fl^,i] = {e,a^,i, d^i} = D2) while for A^ = 6 
[d^,i] = {e, d^, d^,i, d^i, d'^i} = Z?3. Group -D3, on the other hand, is non-Ahelian (id^ = a^i), unlike group 
D2 (id^ = a^i) and, as a consequence, possesses not only one-dimensional, but also two-dimensional ir- 
reps. This is precisely why w e are not permitted to completely decouple the above linearized systeno 
(see [Chechin fc Zhukovl . l2006l ]^. 

In spite of this difficul ty, we can still simp l ify the linearized system 5 = J{t)5 considerably using group 



theoretical methods (see [Chechin &: Zhukovl . l2006l : IChechin et al. 
of this splitting for the case N = 6: 



20071 1). Let us present the final result 
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-4yi 



0, 



(66a) 
(66b) 



m + 2y3 
2/5 + 2y5 

y4 + 2y4 
2/6 + 2y6 



P{t)y3, 

P{t)ye, 
P{t)y4. 



(66c) 



(66d) 



Here P{t) = e 3 — AA{t)[l + e 3 ], while P{t) is the complex conjugate of P{t). 

The st ability of the vr-mod e (the bush B[d^,i]) has a lready been extens i vely discussed i n a number 



of papers IBudinskv &: Bountis 



Poggi fc Ruffd . 119971 : lYoshimura 



19831: Sanduskv fc Paed Il994l: iFlacE 11996 : 



2OO4I : IShinoharal . I200I KiJ : iLeo fc L^ 



Chechin et al. 2002 



20071 : iDauxois et al. 



by d iff erent methods and with emphasis on different aspects of stability. In the papers [Chechin et al. 



2OO2I I 

modes, in both FPU-a and FPU 



Chechin et al 



2005 



2005 1 



2005( 1 the stability of all one-dimensional and two-dimensional bushes of vibrational 
j3 models was numerically studied. In particular, this study revealed a 



remarkable property of the FPU-a chain: the stability threshold of the 7r-mode turns out to be one and 
the same when it interacts with any other mode of the chain. By contrast, other nonlinear normal modes 
of the FPU-a (or FPU-/?) c hain exhibit very different stability thresholds when interacting with different 
modes [Chechin et all l2005l | . 



3.4. Bushes of NNMs in monoatomic chains 

To study bushes of NNMs in monoatomic chains, we shall choose the complete set of normal coordinates 
(pk as the basis o f the configuration space. Here, we use the normal coordinates in the form presented 
Poggi fc Ruff 0I . I1997I ]: 



m 



^k 



1 



sm 






n + cos 



2TTk 



-n 



n 



1, 



,N 



A; = 0..Ar-l, 



(67) 



where the subscript k refers to the mode and the subscript n refers to the atom. The vectors ipk {k = 
0,1,2, ... ,N —1) form an orthonormal basis, in which we can expand the set of atomic displacements X{t) 



Actually, this fact can be understood if one takes into account that the two-dimensional irrep is contained two times in the 
decomposition of the natural representation of the considered chain. 



N-l 
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corresponding to a given bush as follows: 

^{t) = Y^ MtWk- (68) 

fc=0 

For example, one obtains the following expressions for the bushes B[a^,i] and B[a^, a^i] [see Eqs. (|35]) . (f37|l ]: 
B[a\i\: X{t) = { xi{t),X2{t),-xi{t),-X2{t) \ xi{t),X2{t),-xi{t),-X2{t) |--- } 

= l^{t)0N/2 + ^{t)03N/4, (69) 

B[a\ a^i] : X{t) = Jl{t)^N/2 + ?(i)'^iV/4- (70) 

Prom the complete basis (|67|) only the vectors 

0N/2 = 4^i-h 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, •••), (71) 



if N/4 = 4^ih -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, •••), (72) 

V jV 



(^3;v/4 = ^(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, ••■)• (73) 



contribute to the two-dimensional bushes ()69p . ()70p . which are equivalent to each other and constitute 
an examples of dynamical domains. For the bush B[d , i], we can find the following relations between the 
old dynamical variables xi{t), X2{t) (corresponding to the configuration space) and the new dynamical 
variables n{t), v{t) (corresponding to the modal space): 

^it) = -^[x^{t)-X2{t)l 
y(t) = -^[X,{t) + X2{t)]. 

Thus, each of the above bushes consists of two modes. One of these modes is the root mode {(p^N/A fo'^ 
the bush B[a^ , i] and jPnia for the bush B[a^ g^i] ) , while the other m o de 0n i2 is the secondary mode. Indeed, 
according to [Sakhnenko &: Chechinl . 119931 : IChechin &: Sakhnenkd . Il998i | the symmetry of the secondary 



modes must be higher or equal to the symmetry of the root mode. In our case, as we deduce from Eq. ()72p . 
the translational symmetry of the mode (pM/i is o^ (acting by this element on (172p produces the same 
displacement pattern), while the translational symmetry of the mode (,?jv/2 is a^, which has twice more 
symmetry elements than that of (^jv/4 [see ([7T]) .([72 |) ]. Note that the full symmetry of the modes (,?jv/4 a-^^d 
(pN/2 is represented by [a'^,a^i] and [a^,i], respectively. 

Let us now return to the Hamiltonian of the FPU system written as follows: 

N N N 

H = T + V=-Y,xl + - Y,{Xn+l - Xnf + - ^(X„+1 - Xnf ■ (75) 

n=l n=l n=l 

Here p = 3, 7 = a for the FPU-a chain, and p = 4, 7 = /3 for the FPU-/3 chain, while T and V are the 
kinetic and potential energies, respectively. We assume again p.b.c. 

Let us consider the set of atomic displacements corresponding to the two-dimensional bush B[a^,i] 

X{t) = {x,y, -y, -X I X, y, -y, -x \ x, y, -y,-x \ • • • }, (76) 

where we rename xi{t) and X2{t) from Eq. (|69p as x{t) and y{t), respectively. Substituting the atomic 
displacements from Eq. ()76p into the Hamiltonian ()75p corresponding to the FPU-a chain we obtain: 

T = ^{x' + y% (77) 

V = —(3x2 _ 2xy + 3y2) + _^ (^.S _^ ^2y _ ^^2 _ ^3^,^ (yg^. 

These expressions are valid for an arbitrary FPU-q chain with N mod 4 = 0. The size of the extended 
primitive cell (EPC) for the vibrational state (f76ll is equal to 4a and, therefore, when calculating the 
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energies T and V, we may restrict ourselves to summing over only one EPC. In the present case, Lagrange's 
equations ([15|) can be written as follows: 

' x + {3x-y) + a(3x2 + 2xy - y^) = 0, 



y + {2,y-x) + a{x^ - 2xy - 3y^) = 0. 



(79) 



Let us emphasize that these equations do not depend on the number A^ of the particles in the chain (only 
N mod 4 = must hold). 

Eqs. (179p are written in terms of the atomic displacements x{t) and y{t). From them, it is easy to obtain 
the dynamical equations for the bush in terms of the normal modes fi{t) and i^(t). Using the relations (I74p 
between the old and new variables, we find the following equations for the bush B[a^,i] in the modal space 



4a n 
/i + 4|U 7=z^ = 0, 



D + 2iy 



N 
8a 



fiv = 0. 



(80) 
(81) 



The Hamiltonian for the bush B[a^,i], considered as a two-dimensional dynamical system, can be written 
in the modal space as follows: 



H[a\i] = U[? + z>2) + (2^2 ^ ^2) _ ^^^2_ 



N 



(82) 



Note that dynamical equations for all one-dimensiona l and two-dimensional bushes of vibrational 
modes for the FPU-chains were presented in Chechin et all |2005| | . 



3.5. Stability of bushes of vibrational modes in nonlinear chains 

3.5.1. Consideration in modal space 

Since the term "stability" is often used in the literature with different meanings, let us explain in what 
sense we use it in the present paper. 



The stability of bushe s of modes was discussed, in general, in [Sakhnenko fc Chechinl. 11993 



Chechin &: Sakhnenkd . Il998l |. while in the case of the FPU chains it was considered in Chechin et al\ 
2002, l2005 |. Following these papers, we also discuss here the stability of a given bush of normal modes 



with respect to its interactions |i| with the modes which do not belong to this bush. Let us illustrate this 
idea with an example. 

As was shown above, the two-dimensional bush B[a'^,i] for the FPU-a chain is described by 
Eqs. ([80]) ■ (f8T]) . These equations admit a special solution of the form 

/i(t) / 0, v{t) = 0. (83) 

which can be excited by imposing the initial conditions: /^(to) = /^o 7^ 0, /i(to) = 0, z^(to) = 0, z>(to) = 0). 
Substitution of (183)1 into (jHOl) produces the dynamical equation of the one-dimensional bush B [a? , i] (see 



Chechin et a/.l . 120051 ]) consisting of only one mode ^{t) 

i2 + Afi = 



(84) 



with the simple solution 

/i(i) = /iocos(2t). (85) 

taking (with no loss of generality) the initial phase to be equal to zero. Substituting ([85]) into Eq. ([8T]) . we 



obtain 



^ + 



SafiQ 



cos(2i) 



0. 



(86) 



There is an essential difference between the interactions of the modes that belong and those that do not belong to a 
give n bush: we speak about "fo rce interaction" in the former case and about "parametric interaction" in the last case 
(see [Chechin fc Sakhnenkd [T99i |). 
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This equation can be ea sily transformed into the standard form of the Mathieu equation 
Abramowitz fc Stegunl . Il965l | 



i>+[a-2qcos{2t)]L' = 0. (87) 

As is well-known, there exis t domains of stable and unsta ble motion of the Mathieu equation (I87p in 



the a — q plane of its parameters [Abramowitz Sz Stegunl . Il965l | . The one-dimensional bush B [a , i] is stable 



for sufficiently small amplitudes no of the mode //(t), but becomes unstable when ^o > is increased . This 
phenomenon, similar to the well-known parametric resonance, takes place at /^q values which lie within the 
domains of unstable motion of Mathieu's equation (186p . 

The loss of stability of the dynamical regime ()83p . representing the bush B[a^, i]), manifests itself in the 
exponential growth of the mode i^{t), which was identically zero for the vibrational state (j83p and oscillatory 
for the stable regime of (I87p . As a result of z^(t) ^ 0, the dimension of the original one-dimensional bush 
B[a^,i] increases and the bush is transformed into the two-dimensional bush B[a^,i]. This is accompanied 
by the breaking of the symmetry of the vibrational state (the symmetry of the bush B[a^, i] is twice as high 
as that of the bush B[d^,i]). 

In general, we may view a given bush as a stable dynamical object if the complete collection of its 
modes (and, therefore, its dimension) does not change in time. All other modes of the system-according 
to the definition of a bush as a full collection of active modes-possess zero amplitudes and are therefore 
called "sleeping" modes. If we increase the inte nsity of bush vibratiori s , sorn e sleeping modes (because of 
parametric interactions with the active modes Chechin &: Sakhnenkd . 119981 ]) can lose their stability and 



become excited. In this situation, we speak of the loss of stability of the original bush, since the dimension 
of the vibrational state (the number of active modes) becomes larger, while its symmetry becomes lower. 
Thus, as a consequence of stability loss, the original bush transforms into another bush of higher dimension. 
In Section IV we will discuss in more detail the problem of stability of invariant tori from a more general 
perspective, when we study quasiperiodic orbits in which all modes may be active. 

Note, however, that what we described above is the loss of stability of the (one-dimensional) bush 
B[a^,i] with respect to its transformation into the (two-dimensional) bush B[a^,i], using Floquet theory 
and analyzing a simple Mathieu's equation. In the case of a more general perturbation of the A^ particles 
though, one must examine the stability of the bush B[a^,i] with respect to all other modes, as well. This 
cannot be done using Floquet theory as the variations that need to be studied are not about a periodic but a 
quasiperiodic orbit with two rationally independent frequencies. The stability analysis of a multidimensional 
bush is described theoretically in the next subsection and is performed efficiently and accurately employing 
the GALI numerical method outlined in section IV. 

Indeed, as we shall show in Section IV, it is possible that tori may continue to exist at energies higher 
than the thresholds (obtained by Floquet theory) where their "parent" modes become linearly unstable. 
This ensures stability over a wider domain than just an infinitesimal neighborhood of the simple periodic 
orbits (one-dimensional bushes) representing the NNM's of the FPU chain. 

3.6. Stability of bushes of NNMs in the general case 

3.6.1. Some general comments 

Some aspects of bush stability analysis have already been discussed in the previous sections. In this last 
subsection of Section III, we would like to complete the discussion of bush stability using the apparatus of 
the irre ducible representat i ons o f symmetry groups. 

In [Chechin &: Zhukovl . l2006l | . a general group theoretical procedure was developed for simplifying the 



linear stability analysis of periodic and quasipereodic nonlinear regimes in A'^-particle mechanical systems 
with arbitrary groups of discrete symmetry. 

This procedure allows us to split the linearized equations (near a given dynamical regime) into a 
number of independent subsystems whose dimensions can be much smaller than that of the full system. 
The basis vectors of the irreducible representations of the appropriate groups are needed to perform this 
splitting explicitly. On the other hand, we can obtain more easily very useful "splitting schemes", which 
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determine how many subsystems of different dimensions can appear as a result of the above sphtting with 
the aid of the group-representation characters only. Let us discuss these problems in more detail. 

In general, the study of stability of periodic and, especially, quasiperiodic dynamical regimes of me- 
chanical systems with many degrees of freedom presents considerable difficulties. Indeed, for this purpose, 
we often need to integrate large linearized systems of differential equations with time-dependent coeffi- 
cients. In the case of a periodic regime, one can use the Floquet approach, which requires integration over 
only one time-period to construct the monodromy matrix. However, for quasiperiodic regimes this method 
is not applicable, and one needs to solve a system of several differential equations for very long times in 
order to detect stability (especially, near an instability threshold). 

In such a situation, a decomposition (splitting) of the full linearized system into a number of indepen- 
dent subsystems of small dimensions proves to be very useful. Moreover, this decomposition can provide 
valuable information regarding the specific degrees of freedom, which are responsible for the first desta- 
bilization of a dynamical regime as the energy (or a parameter value) is increased. We remark that the 
number of such "critical" degrees of freedom can frequently be rather small. 

3.6.2. Main theorem of bush stability 

Let us consider an A^-degrees-of- freedom mechanical system described by the set of autonomous differential 
equations 

X = F{X), (88) 

where the configuration vector X = {xi{t),X2{t), . . . , X]\f{t)} determines the deviation from the equilibrium 
state X = {0, 0, ... , 0}, while the vector-function F{X) = {/i(X), f2{X), . . . , f]\f{X)} provides the right- 
hand-side of the dynamical equations. 

We assume that Eq. (I88p is invariant under the action of a discrete symmetry parent group Gq. This 
means that for all g ^ Gq Eq. (I88p is invariant under the transformation of variables 

X = gX, (89) 

where g is the operator associated with the symmetry element g of the group Go by the conventional 
definition. 

Using dlH]) and ([89]), we may write X = g'^X, g'^X= F{g-'^X), and finally 

X= gF{g-'X). (90) 

On the other hand, renaming X from Eq. (I88p as X, one can write X= F(X). Comparing this equation 
with Eq. (HOD, yields F{X) = gF{g-^X), or 

F{gX) = gF{X). (91) 

This is the invariance condition of the dynamical equations (I88p under the action of the operator g. Note 
that it must hold for all g G Gq, hence it is sufficient to consider such equivalence only for the generators 
of the group Go . 

Let X(t) = C(t) be an ?TT--dimensional specific dynamical regime in the considered mechanical system 
that corresponds to the bush B[G] (G C Go). This means that there exist functional relations between 
the individual displacements Xi{t) [i = 1,2, ...,A^), which reduce system (j88p to m ODEs in terms of 
independent functions denoted by a{t), b{t), c{t), etc. in the previous section. 

The vector G(t) is, therefore, a general solution of the equation [see, Eq. (JA.ISP in Appendix A] 

GX = X, 
where G is the symmetry group of the given bush B[G] (G C Go). 
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Suppose now we wish to study the stability of the dynamical regime C(t), corresponding to the bush 
B[G]. To this end, we must hnearize the dynamical equations (188p in a vicinity of the given bush, or more 
precisely, in a vicinity of the vector C{t). Thus, let 



X = Cit) + 5{t), 



(92) 



where 6{t) = {6i{t), . . . ,6N{t)} is an infinitesimal A^-dimensional vector. Substituting X{t) from ([92|l 
into ([881) and linearizing these equations with respect to 6{t), we obtain 



6 = J[Cit)]6, 
where ■y[C(t)] is the Jacobi matrix of the system ([55|) : 

dfi 



(93) 



J[Cit)] 



dxj 



X=C{i) 



In (Chechin &: Zhukovl . |2006| | , we have proved the following theorem: 

Theorem 1. The matrix J'[C{t)] of the linearized dynamical equations near a given hush B[G\, determined 
by the configuration vector C{t), commutes with all matrices M{g) (g £ G) of the natural representation 
of the symmetry group G of the considered bush: 

M{g)J[C{t)]=J[C{t)]M{g). 



We introduce, hereafter, a simpler notation for the Jacobi matrix: 



J 



G{t) ^Jit). 



(94) 



Remark. Theorem 1 suggests that if we take a symmetry element g £ Gq that is not contained in G 
{g £ Go\G), the matrix Aig corresponding to g may not commute with J{t). 

Based on the above, we can now apply the well-known Wigner theorem (see, for example, 

Elliott &: Dawbed . |l979t]) to split the linearized system S = J{t)5 into a number of independent sub- 



systems. In fact, it is useful to formulate this theorem in a way that is convenient for our present purposes. 
Consider a reducible representation F of the group G which can be decomposed into a direct sum of 
the irreducible representations (irreps) Vj of this group: 



^^rnjVj 



(95) 



Here mj is the number of times that Tj enters into this decomposition (subduction frequency). We 
denote the dimension of the irrep Tj by Uj . Then the Wigner theorem asserts the following: 

Theorem 2. Any matrix % commuting with all the matrices of a representation T of the group G can be 
reduced to the block- diagonal form 



n 



E"^. 



(96) 



SO that: (a) the dimension of the each block Dj is equal to rujUj and (b) the block Dj consists of sub-blocks 
representing matrices proportional to the identity matrix X„ of dimension rij which are repeated rrij times 
along the rows and columns of the block Dj . 



The structure of one such block Dj = D characterized by the numbers n,- 



D 



I fJ-u'^n fJ-l2^n ■ ■ ■ fJ'lm'^n \ 
\ l^ml-'-n l^m2-'-n ■ ■ ■ l^mm-'-n J 



m is as follows: 



(97) 
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Let us apply the Wigner theorem for sphtting the Unearized system 5 = J{t)5 near the dynamical 
regime (a bush of modes) with symmetry group G. To this end, we will assume that % and F of the theorem 
are respectively the Jacobi matrix J{t) and the natural representation of the group G. 

To implement this splitting explicitly we must pass from the old basis ^oid = {ei, 6*2, . . . , e/v} of the 
natural space to the new basis ^new = {'^i; ^2^ ■ ■ ■ ■, 'Pn} formed by the complete set of the basis vectors 
(fik {k = 1,2, . . . , N) of all the irreps of the group G. If ^new = S^oidj then the unitary transformation 

Jnewit) = S^JoldS (98) 

produces the above-discussed block-diagonal matrix Jnewit) of the linearized system 5new = Jnew{t)5new 
(here 5oid = S5„ 



-'new ) 



3.6.3. Splitting schemes 

The above-discussed implicit decomposition of the Jacobi matrix J{t) is a cumbersome procedure and, 
therefore, it is interesting to know beforehand to what extent it can be useful in specific cases. This can 
be easily determined by means of the theory of characters of group representations. 

Let us consider such an approach by determining the splitting scheme of the linearized system 5 = 
J{t)5. In other words, we would like to find out how many subsystems of different dimensions can be 
obtained as a result of the decomposition of the Jacobi matrix. 

Each block T)j in the decomposition (I96p of the matrix % = J{t) generates an independent subsystem 
of fijirij equations of the linearized system. However, each subsystem automatically splits into rij new 
subsystems consisting of ruj differential equations, as a consequence of the specific structure of these 
blocks [see Eq. (f97|) ]. For example, if a certain block for the matrix J{t) possesses the form {jij = 3, 
rrij = 2) 



V /^2i2^3 /^222^; 



it is easy to check that one can obtain from the system 5 = J{t)5 three independent 2x2 subsystems with 
one and the same matrix 

_^ ^ / /^ii M12 



Ai21 1^22 

Thus, to obtain the splitting scheme of the linearized system 5 = Jb, for each irrep Fj of the bush 
symmetry group G, one must find two numbers — Uj (the dimension of Fj) and rrij (the subduction frequency 
from Eq. (l95|) ). 

The multiplicity numbers rrij can be found by means of Eq. ()A.13P of Appendix A in terms of the 
characters Xj id) of ths irreps Tj , which are well-known for many groups of discrete symmetry (in particular 
the point symmetry groups) and XT id), which represents the character of the reducible representation 
F. These characters can be determined without the explicit construction of the matrices of the natural 
representation F. 

A quick and simple met hod for obtaining xvid) can be found in many textbooks (see, for exam- 
ple, [Landau Sz Lifshita . ll977l |) in connection with studying the small vibrations of multiatomic molecules. 



The main point of the method is that, for a fixed g £ G, the only atoms that contribute to xr{g) are 
those whose position does not change under the action of g on the molecule. Moreover, this contribution 
is determined by the trace of the three-dimensional matrix associated with the symmetry element g. For 
example, the contribution to xrig) from any atom invariant under the action of a rotation by an angle ip 
is 1 + 2cos(99), while that for a mirror rotation is —1 + 2cos(99). 

Let us illustrate this in the case of an A^-particle monoatomic chain for A^ ^ 1. We denote by B[a"^, . . . ] 
the bush B[G], with the symmetry grou p G containing the transn ational subgroup [a"*] and dots standing 



for other generators of the group G. In [Chechin Sz Zhukovl . l2006l | . the following theorem was proved: 
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Theorem 3. The linear stability analysis of any bush B[a^'^, . . .] in an N -degrees- of- freedom monoatomic 
chain can be reduced to the stability analysis of isolated subsystems of the second order differential equations 
with time- dependent coefficients whose dimension does not exceed the integer number m. 

Corollary. If the bush dimension is d, one can perform the above analysis for subsystems of autonomous 
differential equations with dimensions not exceeding (m + d) . 

For example, the scheme for the linearized system about a one-dimensional bush B[a'^, i] of an arbitrary 
nonlinear chain with A^ = 12 particles allows us to split the 12 linearized equations into two one-dimensional, 
two two-dimensional and two three-dimensional subsystems, which are all independent from each other. 



3.6.4. Stability diagrams for NNMs in the FPU chains 

The above discussed decomposition of the full linearized system 6 = J{t)5 into independent subsystems of 
small dimension permits one to analyze efficiently the stability of a given bush in monoatomic chains with 
arbitrarily large number of particles. Using this idea, the stability diagrams for all the one-dimensional 
bushes (NNMs) in both FPU-a and FPU-/? chains were obtained in Chechin et 0^ . 12005 ]. As an example. 



we reproduce here in Fig. [6] the stability diagram for the bush B[d^,in] of the FPU-/3 chain. 



-Ai 



t 



M 



Fig. 6. Regions of stability (white color) of different modes of the FPU-/3 chain, interacting parametrically with the one- 
dimensional bush B[a , m]. 



In this diagram, each point (A, q) provides a value of the bush mode amplitude A and the wave number 
q = -^ associated with the index j of a particular mode. The black points {A, q) correspond to the case 
where the mode j = ^2^ becomes excited because of its parametric interaction with the mode of the 
bush B[a^,in]. The white color denotes the opposite case: the corresponding mode j is zero initially and 
continues to be zero in spite of its interaction with the considered bush. Such diagrams allow one to study 
sta bility of one-dim ensional bushes not only for finite A^, but also for the case A^ — )• oo (for more details 
see 



Chechin et all . [2005j 1. 



In the particular diagram depicted in Fig. [6l we denote the permissible values of the wave number 
q for N = \2hy the dotted vertical lines {q = 2'irj/12,j = 1, 2, ..., 12. The black color corresponds to the 
regions of unstable motion in the plane {q — A). From Fig. [6] one can see that, for an FPU-/3 chain with 
A^ = 12 particles, the NNM B[a'^,iu] with c = {1, -1, -1, 1 | 1,-1,-1,1 | 1,-1,-1,1} is stable for every 
amplitude A > (at least up to ^ = 20), since the vertical dotted lines do not cross the rabbit- ear 
shaped black (unstable) regions. In fact, the mode B[a^, in] is also stable for chains with A^ = 4 and N = 8 
particles, since the corresponding vertical dotted lines are even more distant from each other than for 
N = 12. 

In the case A^ = 16, however, there exist dotted vertical lines which are tangent to the rabbit-ear 
unstable regions, while for A'^ > 16 (note that A'^ mod 4 = must hold!) these lines begin to cross the 
unstable regions. Therefore, we can conclude that the considered nonlinear normal mode in FPU-/? chains 
with A^ > 16 becomes unstable for amplitude values, which lie within the black regions in Fig. [6l The case 
A^ = 16 represents the "boundary" between stable and unstable behavior of the B[a^,iM] NNM. 

Using Fig. [6l one can verify that the critical amplitude Ac of the nonlinear normal mode B[a^, iu], at 
which stability is lost, decreases with increasing A'^. Thus, one can evaluate numerically the corresponding 
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scalin g law of the function AJ N) in the thermodynamic limit, where N ^ oo and i? — )• oo with E/N 



fixed Chechin &: Ryabovl . l201Cll |. It is also interesting to note that if the amplitude of the NNM B[a ,iu] 



exceeds a certain critical value A* w 1.19 (see Fig. [6]), this NNM becomes stable again! Indeed, for A > A*, 
the horizontal lines of A =constant do not intersect any unstable regions. 

Let us summarize the results reviewed of this section concerning the stability of bushes of NNMs: 

(a) If the amplitude of a given NNM (one-dimensional bush) exceeds a certain limit, this mode can lose 
its stability because of specific interactions with some other NNMs of the considered system. As a result, 
the symmetry of the vibrational state becomes lower and the system passes to the dynamical regime of a 
higher dimensional bush whose dimension n > 1. This means that the periodic motion of the considered 
NNM becomes quasiperiodic with n > 1 basic frequencies, which are rationally independent. 

(b) Of course, one cannot use Floquet theory to study stability of bushes whose dimension n > 1, but 
the group theoretical method presented in subsection III.F can simplify significantly the analysis. Indeed, 
the methods we have described allow us to split the linearized equations of motion near a bush into a 
number of independent subsystems whose dimension can be much smaller than that of the full system. 
Thus, when a given bush loses its stability and transforms into another one of higher dimension, one may 
ask: What new modes turn out to be involved into the corresponding vibrational process? 

The following questions also arise: How can one determine the dimension n of these quasiperiodic orbits 
and the frequencies of the corresponding tori? What are the properties of the dynamics in their vicinity? 
How can the stability (and hence the existence) of these tori be determined accurately and efficiently so as 
to gain a more global picture of the stability of motion of the system? To answer these questions, we need 
to examine the phenomena of energy localization and NNM interaction on low-dimensional tori, employing 
analytical methods as well as the numerical approach of the Generalized Alignment Indices (GALI), as we 
describe in Section IV below. 



4. Energy localization on q-tori and stability of motion in FPU chains 

In sections I-III of this paper, we discussed extensively a group theoretical method by which one can con- 
struct exact solutions of mechanical systems described by A^-degree of freedom Hamiltonians. In particular, 
using as an example the FPU A-particle chains, under periodic boundary conditions, we exploited sym- 
metry groups of their equations of motion to determine all continuations of linear normal modes (termed 
here nonlinear normal modes or NNMs), which constitute the so-called one-dimensional bushes. Next, we 
considered linear combinations of these NNMs, in ways which do not excite other modes and were able to 
construct higher (s > 1) -dimensional bushes, which at low enough energies are quasiperiodic, possessing n 
generally incommensurate frequencies. 

In this section, we complete our review by describing two more recent developments in the study 
of lo w-dimensiona l tori of the FPU system, i n wh i ch one of us and his co-wo rkers have played a crucial 
role Skokos et all 120081 . 120071 : iBountis et all 120081 : IChristodoulidi et aU 120101 ]: The first one starts, as in 



the earlier sections, with the construction of s(> l)-dimensional bushes by linear combinations of NNMs 
of low wave numbers g = 1, 2, 3... and establishes the important property of exponential localization of the 
energies of the Fourier modes belonging to the so-called q-tori , through which the well-kn own paradox of 
FPU recurrences can be better understood Bountis et al.l . 120081 : IChristodoulidi et aLl . l2010l |. It is important 



to emphasize that this approach does not depend on the symmetries of the system or the type of boundary 
conditions. 

The second development concerns a more comprehensive method for studying the stability of quasiperi- 



odic orbits, with out resorting to Flq c 



GALI indicators ISkokos et a.li 2008 



uet a nalysis, but using instead the recently introduced spectrum of 
20071 1 . Indeed, the GALI's asymptotic behavior as t — )• oo (typically 



up to t = 10 — 10 ): (i) determines the dimension s of the torus, if the GALI^ indices tend to a constant for 
fc = 2, 3, ..., s (and decay by algebraic power laws for s < k < 2N) and (b) predicts the instability of a torus 
beyond an energy threshold Eth if, for E > Eth, all GALI^ start, at some time t, to decay exponentially, 
indicating that the orbit we are following is not quasiperiodic but diffuses chaotically away from a torus 
that at energies E > Eth no longer exists. 
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4.1. Energy localization on q-tori and FPU recurrences 

We shall use as our main example the FPU Hamiltonian for a chain of N particles, (|99p . which we rewrite 
here for convenience: 

1^1^ a ^ ^ 

^ = 2 IZ^fc + 2 X^^^'^+i ~ ^fc)^ + 3 ^(^k+i - Xkf + ^ ^{xk+i - Xk)"^ (99) 

fc=l A:=0 A:=0 fc=0 

We shall discuss primarily f.b.c, setting xq = xn+i = 0, since the case of p.b.c. does not present any major 
differences . 

For the FPU /3-model (a = 0) of this mechanical system the normal mode canonical variables {Qq, Pq) 
are introduced by the canonical transformations 



5=1 ^ ^ g=l ^ ' 

Substitution of (jlOOp into (j99p brings the Hamiltonian (|99p in the form H = if 2 + -^4, where the quadratic 
part represents A^ uncoupled harmonic oscillators 

i^2 = E ^ 2^ ^ (101) 

q=l 

with normal mode frequencies 

^9 = 2sin(^^^|^), l<g<iV . (102) 

On the other hand, the quartic part of the Hamiltonian becomes 

q,l,m,n=l 

where the coefficients Cq^i^rn,n are non-zero only for particular combinations of the indices q,l,m,n, i.e. 

„ _Jlif(7ib/ibm,±n = . .^ 

<^q,l,m,n -|_lifg±/±j„±„ = ±2(7V + 1) ^^^ 

with all possible choices of the it signs taken into account. Thus, in the new variables, the equations of 
motion are expressed as follows: 

/3 ^ 

Qq + ^^qQq ^ ~o71\/~j~TT / ^ ^q,l,m,n^'q^h^^m^'nQlQmQn ■ (105) 

i,m,n=l 

The question we wish to pursue here is what happens to the (linear) normal modes 

Eq = {P^ + nlQl)/2, q = l,2,....,N (106) 

when /3 7^ 0? In other words, if we choose initial conditions Xk = Xk{0), Pk = Pk{0), k = 1,2,.. .,N in 
(I loop such that a small number (one, two, or more) modes ()106p are excited, with energies Eg, q = 1,2, .., 
will these modes retain the total energy E (and only exchange it among themselves), or will they share it 
equally with all other modes, as statistical mechanics of ergodic systems would require? 

This was precisely the ques tion asked by Fermi, Pasta and Ulam in their famous experiments of the 
early 1950's Fermi et all Il955l |. As is well-known, they obtained the surprising result that, at small E 
values, energy was not shared among all modes, but kept returning to the initially excited ones, giving rise 
to the so-called FPU r ecurrences. There ha.ve be en hundreds of papers to date analyzing this phenomenon 
(for a recent review se e [Berman &: Izrailev . l2005| ) and we shall also provide our own explanation here, using 
the concept of q-tori Christodoulidi et all \201(J\ . The difference between Q'-tori and the multidimensional 
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bushes discussed in earlier sections is that the set of NNMs excited in a bush does not change in time, 
while in g-tori all modes (|106p are active in gen eral, albeit with e xponentially decreasing energies. 



We now recall the theorem by Lyapunov [Lyapunovl . Il892[ | mentioned in the Introduction, which 
states that if the frequencies O^ of the linear normal modes (jl02p of an A^-degree-of-freedom Hamiltonian 
system are incommensurate, these modes survive for /3 7^ 0. This condition of non-commensurability is 
indeed satisfied ii i the f.b.c . case ()102p for many N, e.g. a power of two minus 1, or a prime number minus 



1, among others [Hemmeij . Il959l |. It fails, however, in the p.b.c. case, for which the frequency spectrum 
given by 

f^g = 2sin(^), l<q<N . (107) 

is commensurate for every N, since sin{q7r/N) = sin{{N — q)TT/N), for q = 1, 2, ..., [A^/2]. Thus, the 
arguments and methods based on discrete symmetry groups, which we presented in section II, turn out to 
be extremely useful for proving the existence of NNMs in the p.b.c. case. 

The group theoretical approach considers orbits as dynamical states, see Eq. ([7]), consisting of pat- 
terns (extended primitive cells) which are periodically repeated, see e.g. p^ - (pO]) . ([30|) . By construction, 
therefore, these states are built by modes that correspond to high wavenumbers q = N/2,3N/4:, ..., etc., 
while their symmetries may be expressed by A^~dimensional vectors of the type shown in (j67p . Indeed, 
these vectors provide an orthonormal basis, by which one can construct multi-dimensional bushes using 
linear combinations of such NNMs, as described in subsection III.D. 

There is, however, a more direct way to construct quasiperiodic orbits, which neither employs sym- 
metries nor uses special g-values to achieve this purpose: It is the well-known Poincare-Linstedt method, 
which starts with any number of modes with incommensurate frequencies ujq., for i = 1,2, ...,n, and ap- 
proximates the orbits through the solutions of the linear problem, Xi{t) = AiCOs{ujq.t). Then, through 
interactions provided by the nonlinear terms of the equations of motion, new q'-modes are born, which, 
interestingly enough, have exponentially decaying energies in the Eg vs. q space. These so-called g-tori are 
seen to be n-dimensional (with n as low as 2 or 3), depending on how many modes are initially excited 
and are stable in the sense that the tangent space of the corresponding quasiperiodic orbits is spanned by 
n linearly independent vectors, at least for very long times and sufficiently low energies. 

Let us see how all this is done explicitly on an example of the FPU-/3 model ()99p with f.b.c. and N = 8: 
We shall construct a 2-dimensional torus representing the continuation of the gi = 1 and q2 = 2 linear 
normal modes, Qi{t) = Ai cosflit, Q2 = A2 cosil.2t, for /3 7^ and special values of Ai,A2. Following the 
Poincare - Lindstedt approach, we seek solutions Qq{t), q = 1, . . . ,8 expanded as series in the parameter 
a = (3/2{N + 1), namely 

Qgit) = Q^^Ht) + aQ^;\t) + a^Q^^\t) + ..., q = l,...8 . (108) 

For the motion to be quasi-periodic on a two-torus, the functions Qq (t) must, at any order r, be trigono- 
metric polynomials involving only two frequencies oji and UJ2 (and their multiples). Moreover, these must 
arise as small corrections of the linear normal mode frequencies Qi, O2) given by series expansions in powers 
of a, as: 

w, = f^g + (To; W + a^wf ) + . . . q = l,2 . (109) 

As is well-known, the ujq , i = 1,2,... in (jlOSp are determined by the requirement that no secular terms 
(of the type tsmujqt etc.) arise in the equations of motion. 

For example, the equation for the ffist mode, contains on its right hand side (r.h.s) terms of the form: 

Qi + ^iQi = -a{2,^\Q\ + ml^lQiQl + Snln^Q^Qs + ...) (110) 

Substituting fli on the l.h.s. of Eq. (|110p by the expression obtained by squaring (jl09p and solving for Qf, 
inserting the expansions (jlOSp into (jllOp and grouping together terms of the same order, we find that at 
zeroth order Q\ -\- cofQi = 0, while at first order the following equation is satisfied 

- snfn^iiQf^fQP + ... (Ill) 
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Repeating the above procedure for modes 2 and 3, we find tliat their zeroth order equations also take the 
harmonic oscillator form, i.e. we have: 

Qf^+cofQf^=0, . = 1,2 gf +OiQf =0 . (112) 

Note that the corrected frequencies uji, 0J2 appear in the zeroth order equations for the modes 1 and 2, while 
the uncorrected frequency Vt^ appears in the zeroth order equation of the mode 3 (similarly, f24, . . . , Jig 
appear in the zeroth order equations of the modes q = 4, ..., 8). We now proceed with the solutions of ()112p 
(with zero velocities at t = 0) which read: 

q[°) (t) = Ai cos wii, , qf' (t) = A2 cos uj2t, Qf^ (t) = A3 cos Qst 

where the amplitudes Ai,A2, A3 are arbitrary. If the orbit is to lie on a two-torus with frequencies wi 
and UJ2 only, we must set A3 = 0. In the same way, the zeroth order equations Qq + fi^Qq = 0, for 
the remaining modes q = 4,..., 8, would yield solutions with new frequencies and hence we need to set 
^4 = ^5 = ... = ^8 = also. Thus, we are left with only two non-zero free amplitudes to determine: Ai, 
A2. 

As we explain in detail in [Christodoulidi et all |2010t], equation (jllip can be simplified dramatically 



and upon substitution of Q\ (t) = Ai coswii, Q2 {t) = A2 cosuj2t reduces to: 



Q\' +LoJQ\' =2niu}\>Aicoscjit-3nfAlcos^ujit 

— 6r2i02^i^2 coscuitcos uj2t . (113) 



We now fix u:\ so that no secular terms appear in the solution, yielding 



4'^ = l^l^l + IaI^i^I 



while we also determine Q\ as: 



(1) _ 3Af Of cos 2,0Jit ^AiAl^l^l cos{uJi + 2uj2)t 
^1 (*) ~ 00 2 '" 



32a;f 2[{uJi + 2u}2f - ujj 

3AiA|Ofr2|cos(a;i - 2u;2)t 
2[{uji - 2002)^ - cof] 



^ o^l^2-.1..2^..ovu.i-.u.,;. ^^^^^ 



Repeating the above analysis for the second mode, Q2{t), we obtain the frequency correction of the second 
mode: 



4^) = ^-Al^l + ^Alnl^2 



and the first order solution 



(1) _ 3A^0^ cos 3uj2t 3A\A2^\^l cos(2a;i + oj2)t 
Q2 (*) - 57^:2 ^ 



32a;| 2[{2uji + UJ2Y - ^l 

3Af Aaflf 0| cos(2a;i - ^2)* 

2[{2u:i - UJ2? - ujI] 

which has a similar structure as the first order solution of the first mode. For the third order term there is 
no frequency correction, and we find 



+ ^77:7- , . ^2_„,21 ill^i 



(1) , , _ Afilfr^s /3cosci;it cos3a;it 



4 \ujI-^1 9w2_f^2 

?,AiA^Q.iVL^Q.3 / cos{(jJi — 2uj2)t cos(a;i -|- 2a;2)i 2cosa;iA 

^ 4 V (^1 - 2W2)2 - ^i ^ (Wl + 2uJ2f - 02 + ul -^l) • ^^^^' 

We thus proceed up to the sixth mode and find that all the functions Q3 , ■ ■ ■ ,Q\, are equal to zero, while 
the functions ^3 , • • • , Qg are non zero. However, a new situation appears when we arrive at the seventh 
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Fig. 7. Comparison of numerical (points) versus analytical (solid line) solutions, using the Poincare - Lindstedt series up 
to order 0{a ), for the temporal evolution of the modes (a) q = 1, (b) g = 3, (c) g = 7, when Ai = 1, A2 = 0.5, and 
N — 8, 13 = 0.1 and (d) the time evolution of the indices GALI2 to GALIg up to a time t = 10 show that the motion lies an 
2-dimensional torus. 



and eighth modes. A careful inspection of the equation for Qy (^) shows that there can be no term on the 

r.h.s. which does not involve some of the functions Q3 , • • • , Qg , according to the selection rules for the 

coefficients of Eq. (|lU4p . Since all these functions are zero, we must set Qy (t) = 0, so as not to introduce a 
third frequency in the solutions, whence the series expansion (|108|) for Q'j{t) necessarily starts with terms 

of orde r at least Oia"^). The same holds true for the equation determining Qg (t). 

In Christodoulidi et al\ . 120101 ] it is shown that the Poincare-Linstedt scheme is consistent, in the 
sense that the series constructed above will never encounter terms with zero denominator, as long as our 
parameters lie in the complement of a Cantor set in the uji,uj2 (or Ai,A2) space. We then discuss the 
question of convergence of the series, which is in general a very difficult issue. However, as we also explain 
below, the fact that the amplitudes squared (or,equivalently the energies) of the higher order terms of 
our series are seen to decrease exponentially with q suggests that our series may indeed be absolutely 
convergent, after all. 

To verify this numerically, let us use the analytical solutions (|108|) at t = as initial conditions 
to integrate the equations of motion and compare the computed Qi{t), i = 1,3,7 with our analytical 
formulas. We also apply the GALI criterion (see subsection IV. B) to show that the solutions lie indeed on 
two-dimensional tori. As shown in Figures [7^-c, the numerically evaluated Qi(i),Q3(t) and Qjit) are in 
excellent agreement with the analytical solution truncated at second order with respect to a, when A^ = 8, 
/3 = 0.1, and Ai = I, A2 = 0.5. 

Concerning Fig. [TJi, let us note that, if an orbit lies on a stable s-dimensional torus, the indicators 
GALI2,. . . ,GALIs, discussed in the next subsection, are nearly constant, oscillating ab out non-zero v alues , 
while the GALI^+j, j = 1,2,...,2N — s decay asymptotically by power laws, as t~^ jSkokos et all 120071 . 
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Fig. 8. The average harmonic energy Eq of the q-th mode as a function of q, after a time T = 10 for (a) the 2-torus of 
FiglZK fl'iid (b) a 4-torus solution (open circles). The stars are Eq values calculated via the analytical representation of the 
solutions Qq{t) by the Poincare -Lindstedt series. The filled circles show a theoretical estimate based on the average energy of 
suitably defined groups of modes (see Ea. H120[) and relevant discussion in the text). 
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This is precisely what we observe in Fig. [TJi. Namely, the GALI2 index stabilizes at a constant value 
GALI2 — 0.1, while all subsequent indices, starting from GALI3 decay following a power law as predicted 
by the theory. Thus, we conclude that the motion lies on a 2-torus, exactly as suggested by the Poincare- 
Lindstedt construction, despite the fact that some excitation was provided initially to all modes. In fact, 
we have carried out many such experiments, starting initially with s > 2 modes and, in every case, the 
analytical solutions were found to be very accurate, while all the GALI indicators showed that they lie on 
s~dimensional tori. 

Let us now evaluate, for such excitations, the average harmonic energy Eq of the q mode and plot it as 
a function of the wavenumber q, after a time interval T = 10^. As we observe in Fig. [8^,b, the numerical 
results shown by the white circles are very close to the analytical ones represented by the stars, for an 
initial excitation with s = 2 and one with s = 4 modes respectively. Furthermore, what is clearly evident 
in these semi-log plots, is that the magnitude of the Eq falls exponentially, though not with one and the 
same exponent. Let us examine more closely this phenomenon of energy localization: 

The average size of the oscillation amplitudes of each mode on an s-dimensional q-torus follows from 
norm estimates of all terms of the form Q \l) (t) appearing in equations such as Eq. (jll4p . where q^^' denotes 

the indices of all modes excited at order k of the perturbation. If we define by A^'^' the mean value of all 
the norms ||Q„(fc)||, we obtain the estimate: 



Aik) 



{CsfA, 



k A'ik+l 







2k + 1 



(117) 



where Aq = A^^' is the mean amplitude of the oscillations of all modes excite d at the zeroth order of th e 



theory, and C is a constant of order 0(1). The proof of (|117p is provided in Christodoulidi et all l2010l | 



Appendix B, in which the anal ytical estimate C ~ 3/2 was given. In fact, (|117|) is a generalization of the 
estimate given by Flach et al. Flach et all l2005l . l2006l | for q-breathers, while the two estimates become 
identical (except for the precise value of C) if one sets s = 1 in Eq. (|117p . and qg = 1 in Flach's q-breather 
formula. 

Note, of course, that, with |CsAq| < 1, Eq. (J117p clearly implies that the amplitudes (and corresponding 
energies) of the modes of these s-dimensional tori are decreasing ex ponentially with fc . This is reminiscent 
of the localization properties of q-breathers, studied by Flach et al. Flach et aU l2005l . l2006l | in connection 
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with the paradox of FPU recurrences. We, therefore, decided to cah ah s-dimensional tori whose modes 
are exponentiahy locahzed q-tori. 

In fact, these s-dimensional g-tori, when excited at high enough energies, give rise to the so-called 
"natural packets" of ttt. > s modes, which c ause the breakdown of FPU recurrences by becoming resonant, 
as suggested in a number of papers (see e.g. Berchialla et a?.l . |2004| | . Luca et aUll999t |) on the FPU /3-chain. 



Thus, (7-tori offer a "bridge" between the scenario described in these references and the arguments based 
on energy localization of q-breathers to provide a more complete interpretation of the FPU recurrences. 

Furthermore, using ()117p . it is possible to obtain 'piecewise' estimates of the energy of each group, 
using a formula for the average harmonic energies E^ > of the modes q^ ' . To see how this is done, note that 
the total energy E given to the system can be estimated as the sum of the energies of the modes 1, . . . , s 
(the remaining modes yield only small corrections to the total energy), as 

On the other hand, the energy of each mode q^^' can be estimated from 

I02 f /5 ^'^.(.)^2 -'s^cs^r'Ar' 



E(f^) ~ lf^2 ( (^ ) (aW)2 
2 1^^ \2{N + 1)J ^ ' 



22fc+l('j\r_|_ l\2k+2 

which, in terms of the total energy E^ yields 



^,„^E/CW + l)-g) . (jjs) 



s \ vr^s^ 



Eq. (jllSp is very similar to the corresponding equation for q-breathers Flach et all 120061 ] 



ii;(2fc+i),o - E,, ^ ^^-^-3 j (119) 

where go is the unique mode excited at zeroth order of the perturbation theory and the integer s plays in 
Eq. ijllSp a role similar to that of go in Eq. (jll9p . This means that the energy profile of a q-breather with 
qo = s obeys the same exponential law as the energy profile of the s-dimensional q-torus. 

Still, the most important feature of the solutions on g-tori is that their profile remains invariant as N 
increases, provided that: i) a constant fraction M = s/N of the spectrum is initially excited, (i.e. that s 
increases proportionally with N), and ii) the specific energy e = E/N remains constant. Thus, in terms of 
the specific energy e, (IllSp takes the form 

M V 7r4M4 ) ^ ^> 

i.e. the profile becomes independent of N . A similar behavior is recovered in the q- breather solutions 



provided that the 'se ed' mode go varies linearly with A^, as was shown in detail in refs. [Flach et al\ . 12007 



Kanakov et a?.! . l2007l | . 



The 'stepwise' profiles predicted by Eq. (jll8p in the case of exact q-tori solutions are shown by filled 
circles in Fig. [HI concerning examples of a 2-torus (see Fig. [7^) and a 4-torus solution (in all fittings 
we set C = 1 for simplicity). From these we see that the theoretical 'piecewise' profiles yield nearly the 
same "plateaus" and average exponential slope of the profiles obtained numerically, or analytically by the 
Poincare - Lindstedt method. 

Naturally, regarding the relevance of the q-tori solutions for the interpretation of FPU recurrences, one 
must verify whether Eqs. (lll8p or ()120p retain their predictive power in the case of generic FPU trajectories 
which, by definition, are started close to, but not exactly on a q-torus. 

To answer this question let us examine the results presented in Figures [9] and [TOJ Figure [9] shows the 
energy localization profile in numerical experiments in which ^ is kept fixed (/3 = 0.3), while N takes the 
values N = 64, N = 128 and N = 256. In all six panels of Fig. [9] the FPU-trajectories are computed 
starting with initial conditions in which only the s = 4 (for A^ = 64), s = 8 (for A^ = 128) and s = 16 
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Fig. 9. The average harmonic energy Eq of the q-th mode over a time span T — W as a function of q in various examples 
of FPU-trajectories, for /3 — 0.3, in which the s(= N/16) first modes are only excited initially via Qq{0) = Aq, Qq{0) = 0, 
q = 1, . . . ,s, with the Aq selected so that the total energy is equal to the value E = H indicated in each panel. We thus have 
(a) iV = 64, S = 10"", (b) iV = 128, E = 2x 10""*, (c) iV = 256, E ^ Ax 10""*, (d) iV = 64, S = 10"^, (e) iV = 128, 
E = 2 X 10"'^, (f) A'' = 256, E = 4 x 10"'^. The specific energy is constant in each of the two rows, i.e. e = 1.5625 x 10~^ in 
the top row and e = 1.5625 x 10~ in the bottom row. The dashed lines represent the average exponential profile Eq obtained 
theoretically by the hypothesis that the depicted FPU trajectories lie close to q-tori governed by the profile (|120|) . 



(for A^ = 256) first modes are excited at t = 0, with the excitation amphtudes being compatible with the 
values of the total energy E indicated in each panel, specific energy e = 1.5625 x 10~® in the top row and 
£ = 1.5625 X 10~^ in the bottom row of Fig. [9l Thus, these FPU trajectories can be considered as lying in 
the neighborhood of the q-torus solutions, at least initially. The numerical evidence is that if E is small, 
they remain close to the q-tori even after relatively long times, e.g. t = 10^. 

This is evident in Fig. [9l in which one sees that the average energy profiles of the FPU-trajectories 
(at t = 10^) exhibit the same behavior as predicted by Eq. ljllSp . for an exact q-torus solution with the 
same total energy as the FPU trajectory in each panel. For example, based on the values of their average 
harmonic energy, the modes in Fig. [9^ (in which s = 4) are clearly separated in groups, (1 to 4), (5 to 
12) and (13 to 20), etc., as predicted theoretically for an exact 4-torus solution. The energies of the modes 
in each group have a sigmoid variation around a level value characteristic of the group, which is nearly 
the value predicted by Eq. ()118p . Also, if we superpose the numerical data of the three top (or bottom) 
panels we find that the average exponential slope is nearly identical in all panels of each row, as implied 
by Eq. (jl20p . according to which, for a given fraction M of initially excited modes, this slope depends on 
the specific energy only, and is independent of A^. 

Of course, when the energy is increased, the FPU-trajectories resulting from s initially excited modes 
start to deviate from their associated tori, when these become unstable. As a consequence, the energy 
profiles of the FPU-trajectories also start to deviate from the energy profiles of the exact s-tori. This is 
shown in Fig. [TOh by the fact that the profiles of the FPU-trajectories become smoother and the groups of 
modes less distinct, when the energy increases by a factor 50 (with respect to Fig. Ell), for the same values 
of N and /3. We also observe in Fig. [TOk-f the formation of the so-called 'tail', i.e. an overall rise of the 
localization profile at the high-frequency part of the spectrum. This is a precursor of the evolution of the 
system towards equipartition, which ultimately sets in at high enough energy. 

This exponential localization of energy profiles also occurs in FPU particle chains, when p.b.c. are 
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Fig. 10. Same as in FigUJi but for larger energies, namely (a) E = 0.05, (b) E = 0.1, (c) E = 0.2, (d) E = 0.3, (e) E = 0.4, 
(f) E = 0.5. Beyond the threshold E ~ 0.05, theoretical profiles of the form (|118|l yield the correct exponential slope if s is 
gradually increased from s = 4 in (a) and (b) to s = 6 in (c), s = 7 in (d), s = 10 in (e), and s = 12 in (f). 

imposed. In fact, the profiles in this case become perfectly symmetric with respect to the middle modes 
if one respects the symmetry of the (linear) normal mode spectrum, (jl07p . and excites equally the q 
mode and N — q mode, which have the same frequency, see Fig. [TTJ The existence and stability of the 
corresponding g-tori at low energies is also evident here and is expected to play an equally important role 
in the phenomenon of FPU recurrences, as in the case of the fixed boundary conditions. 



4.2. Stability of the motion on q-tori and the GALI method 

As we have frequently emphasized, one of the most important properties concerning the dynamics of 
conservative mechanical systems, is whether the motion is "regular" or "predictable" for a given set of 
parameters and initial conditions. Since stable periodic orbits are a set of measure zero, this typically 
implies that the solutions we seek are quasiperiodic in time. In A^-degree-of -freedom Hamiltonian models, 
quasiperiodic orbits are generically located around stable periodic orbits and lie on A^-dimensional tori. 

In particular, when studying the phenomenon of FPU recurrences in subsection IV. A, we realized that 
it is crucial to follow NNMs with low wave numbers q = 1,2,3,.... These are simple periodic solutions 
called g-breathers (due to their exponential localization in g-space), whose initial conditions are very close 
to those of the recurrent FPU trajectories. Owing to this proximity, one might expect that as long as 
g-breathers are stable recurrences will continue for arbitrarily long times. 

In fact, the situation is somewhat more complicated. The reason is that exponential localization in 
g-space is also a property shared by the tori surrounding these NNMs, which we called g~tori in subsection 
IV. A. Actually it is their stability that matters and not that of the g-breathers, since, as the energy 
increases, q-tou are seen to persist even after the NNMs at their midst have turned unstable! 



Of course, the stability of q-breathers can be studied by Floquet theory, as shown e.g. in [Flach et al. 

200a], where it was demonstrated that q-breathers with low q = q^ = 1,2, 3, ... are linearly stable as long 
as 



-Ego < 



vr 



6/3(iV + 1) 



0(1/A^2) 



(121) 
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Fig. 11. The average harmonic energy Eq of the q-th mode as a function of q for the FPU-q model with a — 0.3 and P — 
in (|99p evaluated after a time T = 10 . (a) A q-breather with TV = 16 is shown at energy E = .05 is (b) and a 2-torus and 
4-torus solution with A'^ = 31 are depicted at energy E = .001 . Note that in (a), when only the q ~ 1 mode is excited, the 
profile is not as symmetric as it becomes when both q — 1 and q = N — 1 modes are excited equally. 



This result is obtained by analyzing the eigenvalues of the monodromy matrix of the linearized equations 
about a q-breather constructed by the Poincare - Lindstedt series. Furthermore, examining also the asso- 
ciated eigenvectors one gains a complete understanding of the tangent dynamics in the vicinity of these 
simple periodic orbits. 

In the case of q-tori, however, the above techniques are no longer available. Stability is not decided by 
the eigenvalues of a solution matrix and one must seek different ways of exploring the variational equations 
linearized about quasiperiodic orbits. One analytical approach for example, described in section III, is to 
use discrete symmetries of the system to write the linearized equations in the simplest form, expanding 
the solutions in terms of linear combinations of NNMs, so that they "uncouple" as much as possible. 

Interesting as this approach may be, a more accurate stability analysis should focus on the full set of 
variational equations and examine the behavior of their solutions for arbitrary choices of initial conditions. 
This can be achieved by employing t he recently p r opose d method of the Generalized Alignment Indices 



200i]): 



(GALI), which we outline below (see jSkokos et al. 

Consider an orbit 7(t) = (p(t),q(t)), representing a solution of the Hamiltonian equations of motion 
([1]) with initial conditions (p(0),q(0)). The variational equations along this orbit describe the dynamics 
on the tangent space of the solutions of ([1]) and are given by the vector field 



^ = J-M(p(t),q(t))w 



(122) 



where J 



O In 
-In O 



In the N x N unit matrix and M(p(t), q(t)) is the Hessian matrix of the Hamil- 



tonian function evaluated along 7(t), i.e. 

M(7(t)) 



5p9q '^^ ' 



(123) 



Let us now solve equations (11221) for ■wi(i), ■W2(t), ..., Wfc(t), using as initial conditions k randomly 
chosen orthocanonical vectors and construct the unit vectors 



Wi 



W,; 



W,; 



l,...,k 



(124) 
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The Generalized Alignment Index GALIfc(t) is defined at time t as the volume of the parallelepiped pro- 
duced by these k unit deviation vectors Wi, i = 1, ..., k expressed, at time t, by the wedge product 

GALIkit) =\\wiAw2r\ ■■■ AwkW (125) 

As is well-known, this wedge product can be written with respect to the basis ej^. in the form 



Wl A'W2 /\ ■ ■ ■ /\Wk 



Win wii2 

W2ii W2i2 



Wlik 
W2i, 



Cji A ej2 A • • • A ii^ 



(126) 



E 

l<ii<i2<---<ik<2N 

Wkii Wki2 ■ ■ ■ Wki^ 

From this it is easy to see that if at least two of the normalized deviation vectors Wi, i = 1,2, ... ,k are 
linearly dependent, all the k x k determinants will become zero and the volume will vanish. Clearly, the 
norm of this quantity, ||i()i(t) A W2{t) A • • • A Wk{t)\\, is obtained by evaluating the sums 



GALIkit) 



Wli^ Wli2 
W2i^ W2i2 



W2i, 



1/2 



(127) 



E 

l<ii<i2<---<ik<2N 

Wkii Wki2 ■ ■ ■ Wki^ 

As was described in detail in a number of papers Skokos et aLl . 120071 . 120081 : IChristodoulidi &: Bountisl . l2006t |. 
it is possible to analyze the asymptotic behavior of these sums by estimating the largest determinants as 
t — )• oo, both for chaotic as well as quasiperiodic motion. Thus, the following asymptotic estimates are 
derived: 

I) For chaotic orbits: 

-[(CTi-o-2) + (o-l-cr3)H h(o-i-crfe)]i 



GALIfc(t) oc e 



where 



fl > ^2 > • • • > 0"Af-l > <^N = (^N+1 = > aN+2 > • • • > 0"2Ar. 

are the Lyapunov exponents of the orbit Skokosl . l2010l | . 

II) For quasiperiodic orbits lying on an s-dimensional torus: 



(128) 



(129) 



GALIfc(t) 



constant if 2 < A; < s 
-^ \is<k<2N -s 

1 if 2A^ - s < A; < 2A^ 



(130) 



£2{k~N) 

In the "generic case" s = N , (jl30p implies that the GALI^ remain constant for 2 < k < N and decrease 
to zero as ~ \/t^^^^^^> for N < k < 2N. However, as we have seen in this review, the "non-generic" case 
of low-dimensional tori, with s << A^, is also of great physical interest and occurs quite frequently in 
mechanical systems of the type considered here. 

The accuracy of the above formulas, when compar ed against numerical experiments, is remark- 



able. Indeed, as has been demo nstrated in many papers [Skokos et all 120071 . 120081 : iBountis et all l2008l : 
Christodoulidi &: Bountisl . l2006l | the GALI^ reach their asymptotic values after relatively short times. In 
the chaotic case, the indices follow very closely the theoretical values even when the Lyapunov exponents 
(and their differences) are very small and have not as yet converged to their limiting values. Indeed, it is 
important to note that what makes the GALI method so efficient in distinguishing chaos from order is the 
fundamental difference between exponential and power law decay. In the case of chaotic orbits, all indices 
decrease exponentially and chaos is most quickly detected by inspecting the ones with high k values, while 
for quasiperiodic orbits all GALI^ decrease by power laws, except for the first few which are nearly constant 
and signify the dimension of the torus. 

Concerning the calculation of the GALIfc, it is clear that, for large A^, the computations become 
prohibitively time-consuming, due to the number and size of determinants in (jl27p . There is, however, a 
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much faster way to proceed: Observe that the normahzed deviation vectors Wi, i = 1,2, ... , 2N, discussed 
above, when expressed in terms of the usual orthonormal basis of the n-dimensional Euchdian space can 
be written as 



Wi 



2N 

E 



wae 



tj^j , 



1,2, 



,k 



(131) 



Thus, we can write equations ()13ip in matrix form as: 



Wl 
W2 



Wii U'12 ••• Wi2N 
W21 W22 • • • W22N 



Wkl Wk2 



Wk2N 





ei 




ei 




62 


= A- 


62 




_e2N _ 




_e2N _ 



(132) 



Recah that GALI^ measures the volume of a fc-p arallelepiped Pk haviiig as e dges the k unitary deviation 
vectors Wi, i = 1, . . . ,k. This volume is given by ( [Hubbard fc HubbardLll999l |): 



vol(Pfc) = Vdet(A- A^) 



(133) 



since det(A • A ) is equal to the sum appearing in ()127p by virtue of Lagrange's identity [Bourbakil . Il958l | 
(where ( ) denotes transpose). Thus we obtain a much more efficient way of computing GALI^ by the 
formula: 



GALK 



det(A • A^ 



(134) 



where only the multiplication of two matrices and the square root of one determinant appears. 

Let us see how one can apply the above spectrum of indices to the problems considered in this review. 
First of all, we have already seen the usefulness of the GALI^ in determining the dimensionality of the 
g-tori described in subsection IV. A, see Fig. [7li and the discussion below it. Indeed, the near constancy of 
the GALIfc, for k = l,2,...,s, at least for the time intervals studied, is always seen to characterize stable 
tori formed by exciting s linear normal modes of our FPU chains, at low energies. 

Quasiperiodic tori, however, can b ecome unstable by losin g their smoothness and eventually breaking 
down into sets of points called cantori MacKav fc Meissl . Il987l | . This means that the variational equations 
about them possess a non-empty unstable manifold through which small variations about the torus are 
expected to grow exponentially. It follows, therefore, that beyond the critical energy at which a g-torus 
becomes unstable, all deviation vectors lying on the unstable manifold will tend to align in the direction 
with the largest exponent, making all GALI indices fall exponentially fast. Thus, the instability of the 
torus can be verified, e.g. following the k = 2 index, GALI2, which is the fastest to compute, even though 
the indices for higher k values tend more quickly to zero. 

Using this criterion, it is easy to test the stability of q-tori and determine approximately the value of 
the critical energy Ec at which they become unstable. Fig. [12] shows an example of such a calculation using 
trajectories started close to a 3-torus (upper curve), 2-torus (middle curve), together with one starting 
close to a g-breather (lower curve) for fixed /3. Observing the dependence of this instability threshold on 
N, we discover that the Ec values of at which the g-tori destabilize are significantly higher than those 
corresponding to the destabilization of g-breathers. 

This is also inferred from the fact that the square and triangle data in Fig. [12] show a much weaker 
dependence of ^ on A^ than the N^^ law predicted by (I12ip for q-breathers. The relation A oc N~^ is 
seen to hold well only for FPU-trajectories started close to a q-breather with q = 1 (filled circles), whose 
numerical curve is close to the dashed line predicted by Eq. ()12ip . These results strongly suggest that the 
critical thresholds Ec at which g-tori (obtained by exciting the low modes) break down are significantly 
higher than the energies at which the corresponding q-breathers become unstable. 

We, therefore, conclude this section by emphasizing that perhaps the most important feature concerning 
the dynamics of orbits forming q-breathers or q-tori is the exponential localization of their energies, which 
persists even after these solutions have become linearly unstable. This behavior is exemplified in Figure 
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Fig. 12. Fixing /3 we plot in the upper curve (squares) the dependence of the mean oscillation amplitude A on A'^, for critical 
energies Ec determined by the destabilization of a S-dimensional g-torus for FPU-trajectories identified as quasiperiodic by 
teh GALIjj=const. for k = 1,2,3. The middle curve (triangles) corresponds to a similar calculation for an FPU trajectory 
started by exciting initially a 2-dimensional g-torus, while the lower one (circles) refers to a q-breather solution (go ~ 1). The 
dashed line represents the law of Ea. (|12ip and all calculations are carried out up to i = tmax = 10 . 

[THl where panels (a) and (c) show the time evolution of the GALI2 index for two FPU trajectories started 
in the vicinity of 2-tori of FPU-^ systems with A^ = 32 and N = 128, for /? = 0.1 and tmax = W^ ■ In both 
cases, the energy satisfies E > Ec, as the exponential decay of the GALI2 index evidently shows and the 
corresponding q-tori are unstable. Still, as Figs. [TSb.d clearly demonstrate, the exponential localization of 
the energies, at least for the lower part of the spectrum, holds in this case also. 

This localization property of the g-tori, therefore, provides, in our view, a more complete interpretation 
of the phenomenon of FPU recurrences. If there is still something "paradoxical" about it, it is that it is 
not known to date if it holds to t — >• 00, for small but non-zero energies. What is certainly true (as Figs. 
[T3b.d show) is that the destabilization of the g-tori brings about the delocalization of higher g-modes and 
the break down of recurrences, signifying the eventual equipartition of energy, which is expected to spread 
to all modes at t > Teq{e), where e = E/N is the specific energy and become practically observable when 
e is large enough. 

Appendix A 



5. Elements of group representation theory 

5.1. Natural representation 

We begin with an especially useful representation, called here natural representation of the parent group, 
using as an example the square molecule model, depicted in Fig.[2l We can describe the atomic displacement 
patterns of the vibrational states of this molecule at any fixed time t by eight-dimensional configuration 
vectors of the form 

X = {fi\f2\r3\u} = {a;i,X2|x3,2;4|x5,X6|x7,X8}. (A.l) 

where xi, X2 are displacements of the first atom from the equilibrium position (see Fig. [2]) along the axes 
X and Y, while {x3,Xi}, {x5,xq}, {x7,X8} are those of the atoms 2, 3 and 4, respectively. 

Let us act on an instantaneous vibrational configuration of the molecule by a certain symmetry element 
Qi, see Eqs. (|4T]l . Actually, we must act on the displaced atoms whose coordinates are determined by vectors 
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Fig. 13. (a) Time evolution of the GALI2 index up to i = 10 for an FPU trajectory started by exciting the q ~ 1 and q — 2 
modes in the system A'^ — 32, /3 — 0.1, with total energy E — 2, i.e. higher than Ec ~ 1.6. (b) Instantaneous localization 
profile of the FPU trajectory of (a) at f = 10 . (c) Same as in (a) but for N — 64, E = 1.323 (in this case the critical energy 
is Ec = 1.24). (d) Same as in (b) but for the trajectory of (c). 

Rj = Rj + fj, where R^ = {x^,y^) is the equihbrium position of the j-th atom, while fj = {xj,yj) is the 
local vector which determines the displacement of this atom from its equilibrium position. 

Since giRj = giR^ + girj, one can introduce operators gi, which are induced by the symmetry elements 
u gi £ Go, acting on the eight dimensional configuration vector (jA.ip by the following two step process: 

1) transpose the local displacement vectors fj = {xj,yj) in (jA.ip according to the transposition of 
atom equilibrium positions under the action of the symmetry element g^ ; 

2) act on each local displacement vector fj by the symmetry element g~ . 
For example, we obtain for the symmetry element g2- 



hX = {gif2\gifz\gifi\gifi] = {y2, -X2\y?., -x?,\y4., -Xi\yi, -xi}, 



(A.2) 



Note that symmetry element g acts on the vectors f of three-dimensional Euclidean space, while operator g acts on vectors 
of the eight-dimensional configuration space. In general, the induced operator g, acting on an arbitrary function f{f), is 
determined by the equation gfjf) = f(g~f), as shown in any textbook on the theory of symmetry groups (see, for example. 



aetermmed by tne equation gj(r) = f[g r), ; 
[Hammermeshl . 1 19621 : lElliott fc Dawber[|l979l ] ) 
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taking into account relations §2 = 54 and g4{x,y) = {y,—x) ((72 and 54 are rotations about the Z axis 
through the angles 90° and 270°, respectively). 

The above definition of the operators gi induced by symmetry elements gi allows one to construct 
an eight-dimensional natural representation of the symmetry group Gq = C^v for the considered square 
molecule. To this end, we may choose the "natural" basis <1> = {ei, 6*2, ... , e/v} of orthogonal unit vectors 
Ci (having 1 in their ith entry and zero everywhere else) to describe all possible displacements in the 
configuration space. 

Next, acting by an operator 5 (5 G Gq) on the vector e^, we can write gcj as a certain linear combination 
of all basis vectors: 

N 

gej = Y^ M^j{g)ei, j = l,2,...,N. (A.3) 

i=l 

This equation associates the matrix M{g) = {-Mijig)} with the operator g and, therefore, with the 
symmetry element g G Gq: 

g^g^M{g). (A.4) 

The set of matrices A4{g) corresponding to all g €z Gq forms the natural representation T for our mechanical 
system. 

In the case of the square molecule, A^ = 8, therefore, this representation turns out to be eight- 
dimensional — and represents a collection of 8 x 8 matrices Ai{g) which correspond to all the elements g 
of the parent group Gq . 

Moreover, these matrices satisfy the same multiplication rule as the symmetry elements corresponding 
to them: if gigj = g^^ then M{gi) ■ M{gj) = M[gk)- Among them there is the unit matrix [M{gi) = /], and 
for every matrix there also exists its inverse [M{g~^) = M~^(g)]. Taking into account the associative rule 
for matrix multiplication, we can show that all different matrices belonging to the set {M{g) \\/g G Gq} 
generate a certain matrix group homeomorphic to the parent group Gq (In general, a matrix representation 
can contain a number of identical matrices). 

It is well-known that any finite group possesses an infinite number of representations whose dimensions 
can be arbitrarily large. However, there exists only a finite number of nonequivalent irreducible represen- 
tations (irreps). Let us consider this issue in more detail. 

5.2. Irreducible representations 

Let us first use the basis ^ = {ei, ...,e*/v} to construct the natural representation of the parent group 
Go = 6*4^,. Of course, one may use any other basis of the configuration space for the same purpose. This 
leads to considering a new form of matrices corresponding to symmetry elements g £ Gq as follows: 

MneUa) = S~'M{g)S, V5 e Go. (A.5) 

where 5 is a non-singular matrix. Based on this idea, we may ask: "How can one choose the basis of our 
configuration space to make the matrices M{g) of the natural representation as simple as possible?" 

Indeed, let us recall that if a number of matrices commute with each other they can be transformed to 
the simplest possible — diagonal — form via a transformation of the type (jA.Sh . using a suitably chosen 
matrix S. On the other hand, if the matrices do not commute, they cannot be diagonalized simultaneously 
by transformation of the form ()A.5|1 . i.e. such a matrix S does not exist. 

The natural representation is only one example of a matrix representation of a given symmetry group. 
Continuing our discussion of arbitrary matrix representations of finite groups, we now consider represen- 
tations Tj of the group G by complex n x n matrices Ai{g),g G G. These matrices act on the vectors of 
an n-dimensional complex vector space which is called the carrier space o f L. 

In physics, we usually deal with unitary representations [Hammermeshl . 11962 : lElliott Sz Dawbeii . ll979l |. 



All matrices M.{g) of such representations are unitary matrices, satisfying A4^{g) = M. ^{g),\/g £ G. 
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Definition A.l. Two matrix representations T = { A4{g) \ g £ G } and F' = { A4'{g) \ g £ G } of the 
group G are called unitary equivalent, if there exists a unitary matrix S (i.e. 5^ = S^^), such that all their 
matrices obey the relation 

M'{g) = S^M{g)S, (A.6) 

The transformation (|A.6P appears when we pass from one orthonormal basis of the carrier space of a 
given representation to another. Let us now define the matrix representation T associated with a certain 
basis "i = {0i{f), (p2{^i ■ ■ ■ , 0n{f^} of the carrier space: 

Definition A.2. If 

g$ = M{g)^ (for all g G G), (A.7) 

where operator g acts on the basis vectors <^j(r) as g0j{r) = 0j{g^^r), then 

r = { M{g) I Vg E G } (A.8) 

is a matrix representation of the group G. Note that M{g) denotes the transpose of the matrix M{g). 

Let us also denote by the symbol V the carrier space of the matrix representation T of the group G 
and suppose we can write ^ as a direct sum 

V = f2Vj (A.9) 

i=i 

of a number K > 1 of invariant subspaces Vj, which are independent of each other under the action 
of operators g induced by all the elements of the group G. More precisely the above invariance can be 
explained as follows: Each operator g {g G G) acting on any vector v belonging to Vj transforms it into 
another vector of the same subspace, i.e. gv cannot coincide with a vector of any other subspace in the 
sum ()A.9P different from Vj. 

If we choose an arbitrary basis in the subspace Vj , we obtain a matrix representation Vj of the group 
G which is constructed on this basis according to the general definition ()A.7|A.8p . 

Since the basis of the full space V can be chosen as the union of all subspaces Vj entering Eq. ()A.9P , 
we find that the representation T can be written as a direct sum of all the representations Tj constructed 
using the bases of the individual subspaces Vj : 

N ® 

r = Y. r,- (A.io) 

i=i 
The dimension n of the full carrier space V is the sum of the dimensions Uj of all its invariant subspaces 

N 

n = 'S2nj. (A. 11) 

i=i 

If a subspace Vj cannot be split into invariant subspaces of smaller dimensions, it is called irreducible 
and so is the representation Tj, defined on it. On the other hand, the representation T, defined on the 
full space V is called "reducible", if it can be decomposed into a finite number of K > 1 irreducible 
representations . 

As a consequence of Eq. (|A.10p . a reducible representation T possesses block diagonal form, whose 
blocks turn out to be the individual irreducible representation Tj. For the given matrix Di{g) G Tj, each 

of these blocks, D^Ag), represents an n x n matrix corresponding to g £ G, while i numbers the blocks 
of one and the same dimension (re). Observe that, when a certain irrep Tj is contained in the reducible 
representation T a number of times ruj (the subduction frequency), some blocks may be identical. 
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We may, therefore, rewrite equation (jA.lOp in the form: 



i 



where rrij indicates how many times the irrep Tj is contained in the decomposition of the reducible rep- 
resentation r. Let us note that the basis vectors of different copies of the same irrep Tj are essentially 
different. One should also emphasize that matrices of all irreps Tj must be known to obtain the explicit 
decomposition of the given reducible representation, while for the schemes of the decomposition ()A.12p one 
needs to know only the characters of these irreps. 

Definition A. 3. A character x\^] of the n-dimensional representation is the n-dimensional vector con- 
structed by the traces of all the matrices M (g) of this representation: 

x[T] = ^={TTM{g)\geG}, 

V ll'-'ll 

where ||G|| denotes the order of the group G. 

There exists a powerful group theoretical apparatus using the irreps and their characters for solving a 
wide class of crystal and molecular systems. The characters of the irreps of a given group are orthonormal 
vectors. For calcu lating the subductio i i frequency nij of the ir rep Tj from Eq. ()A.12p . one can use the 
following formula JHammermeshl . Il962l : lElliott &: Dawberi . Il979t | : 



m, = (x[r], ^[r,]) ^jr^Y. xr{9)xAg)- (A.13) 

Here x[r] is the character of the reducible representation F, while xi^j] is the character of the irreducible 
representation Tj {x[^j] is the complex conjugate value with respect to x[rj])- 

5.3. Some properties of irreducible representations and their role in Physics 

Irreducible representations (irreps) possess many r emarkabl e properties wh i ch can be found in any textbook 



irreducible representations (irrepsj possess many r emarRabl e properties wn i cn can be round m any textboo. 
on group theory (see, for example, [Hammermeshl . Il962: El liott &: Dawbed . Il979l |). In particular, it is well 



known that any group G of finite order||G|| possesses a finite number of different nonequivalent irreps, and 
this number is equal to the number of classes of conjugate elements of the given group. The dimension (nj) 
of each irrep (Tj) must be a divisor of the group order ||G||, i.e. ||G|| mod Hj = 0. Dimensions of all the 
irreps of the group G satisfy the Burnside theorem: 

N 

Y^n] = \\G\\. (A.14) 

i=i 
Among irreps of any group G there is the so-called identity representation which associates with every 
element g £ G one and the same one-by-one unit matrix. 

Let us now return to the group G = C^y of the square molecule. The order ||G|| of this group is equal 
to 8 (all its elements are listed in Eq. (l4ip ). According to Eq. (IA.14P there exist only three variants for the 
possible set of dimensions of the irreps of the group G = C4„: 

2^ + 2^ = 8, (A.15) 

1^ -Fl^ + 1^ + 1^ + 1^ + 1^ + 1^ + 1^ = 8, (A.16) 

1^ -h 1^ -Fl^ -Fl^ -F 2^ = 8. (A.17) 

In variant (jA.lSp . there are two 2-dimensional irreps (ni = 2,n2 = 2). This is impossible because any 
group must possess at least one 1-dimensional irrep : the identity irrep. Variant (IA.16P [uj = 1, j = 1...8) 
is also impossible because only for Abelian group all the irreps can be 1-dimensional, while the group 
G = G^y is non-Abelian. Therefore, only variant (IA.17P turns out to be admissible: the group G = G^y 
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possesses five irreps : four 1-dimensional (ni 
presented in Table 1 below. 
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1) and one 2-dimensional {n^ = 2), as 



Table 1. Symmetry elements of the group G and their irreps 
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Let us note that the construction of the irreps for a given group G is not a simple mathematical 
problem. However, the expli cit form of all the irreps of point and space groups can be found in many 
textbooks (see, for example, [K ovalevl . ll993l |). 



The importance of irreducible representations of symmetry groups is illustrated by a famous theorem 
due to Wigner, according to which eigenstates of any quantum system are classified by irreps of the 
symmetry group G of its Hamiltonian. Namely, every energy level of the quantum system is associated 
with a certain irrep Tj of the group G, and the degeneracy of this level is determined by the dimension rij 
of the irrep Tj . 

The normal modes, representing exact vibrational regimes in the harmonic approximation, can be 
classified in the same manner. Indeed, every eigenstate of the force constant matrix /C (see Eq. (I5ip ) 



represented by the natural frequency Uj and the corresponding eigenvector Cj can be associated with a 
certain irrep Tj of the parent symmetry group of the considered mechanical system. In both cases, the 
differential equations describing the physical systems (quantum or classical) are linear. In this regard, it 
is interesting to investigate any possible symmetry - determined classification of exact dynamical regimes 
in nonlinear classical systems with certain parent symmetry groups. This problem is fundamental in the 
theory of bushes of nonlinear normal modes, as we demonstrate below. 



5.4. Selection rules for excitation transfer between modes of different 
symmetry 

Let us return to the discussion of A^-degree-of- freedom mechanical systems with the parent symmetry group 
Go whose vibrational state is described by the A^-dimensional configurational vector X{t). 

As has already been discussed, all the vibrational states can be classified by the subgroups of the 
parent group Gq. Let us consider the vibrational state characterized by a certain subgroup G C Gq. This 
means that the configuration vector X(t) must be invariant under the action of all the operators g €z G: 

gX{t) = X{t), ygeG. (A.18) 

The group G = { g \ \/g £ G } consists of all the operators g which are associated with the elements g of 
the symmetry group G. 

Let us now decompose the configuration space into the minimal subspaces, which are invariant with 
respect to the parent group Gq. Every one of these subspaces Vj is the carrier space of a certain irrep 
Tj of the group Gq. As a result of this transformation, the mechanical representation T turns out to be 
decomposed into a number of irreps: 



r 



r 



(A.19) 



The set of all the basis vectors (p^ of the irreps Tj entering in the direct sum (]A.19j) may be chosen as 
a new basis of the configuration space. Therefore, we can decompose the configuration vector X{t) in the 
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above mentioned basis as follows: 



n 



T 



where Aj(t) is the contribution to X[t) from the irrep F^, while i indicates different basis vectors of the 
dimensional irreducible representation Tj. Here we introduce the "supervector" $-,• = {(p^- ,(p, , . . . ,(p- ^ } 

which determines the set of the basis vectors of the irrep Tj, the vector flj{t) = {^- {t),fi- (t), . . . , ij.- ^ (t)}, 
which determines the time-dependence of the individual vibrational modes, and their formal scalar product 
ijlj{t),^jj. Prom the condition of invariance of the configuration vector X{t) with respect to the group 

G, (jAlS]) . it can be deduced (see [Sakhnenko fc Chechinl . Il993l : IChechin efail . |2007| |) that the following 



invariance conditions also hold regarding the individual irreps Tj-. 

{rjiG)flj=flj. (A.21) 

Here we have introduced the restriction (Tj \, G) of the irrep Tj of the parent group Gq onto the subgroup 
G C Gq. The restriction (Fq i G) is the set of all the matrices of the representation Fq of the group Go 
which are associated only with the elements of its subgroup G C Gq. 

Eq. (lA.2ip provides the important selection rules for excitation transfer between modes of different 
symmetry. Let us consider this point in more detail. This equation indeed holds for every individual 
irrep Tj of the parent group Gq. Thus, from ()A.2ip . we can see that jlj is invariant under all matrices 
Aij{g) G (Tj I G), i.e. flj is the common eigenvector of all the matrices Aij{g), g €z G with unit eigenvalue. 

As the general solution of Eq. (IA.21[1 . the vector flj depends on a number of arbitary constants denoted 
by the letters a,b,c,d, etc., in the examples of section HL To find all modes contributing to a given 
dynamical regime with symmetry group G, expressed in Eq. ()A.20p by the vector X{t), one must solve the 
linear algebraic equations ()A.2ip for each irrep Tj of the group Gq. As a result, the invariant vectors flj 
for some irreps Tj may turn out to be equal to zero and thus not contribute to the considered dynamical 
regime. On the other hand, some nonzero invariant vectors flj for multidimensional irreps can be very 
special due to certain relations between their components (for example, some components are equal to 
each other, or differ only by sign). Naturally, the contributions Aj(t) to X{t) from such irreps Tj also 
possess very special features. 

Actually, Eq. (|A.2ip can be considered as a source of certain selection rules for excitation transfer from 
the "root" mode to a number of other (secondary) modes. Indeed, if a certain mode with the symmetry 
group G is excited initially (called the "root" mode), this group determines the symmetry of the whole 
bush. The condition that the appropriate dynamical regime X(t) must be invariant under the action of the 
above group G leads to Eq. (JA.ISP and then to Eq. ()A.2ip . If the vector flj for a given irrep Tj is the zero 
vector, there are no modes, belonging to this irrep, which contribute to X(t), i.e., the initial excitation 
cannot be transferred from the root mode to the secondary modes associated with the irrep F-,-. 

Note that basis vectors associated with a given irrep Tj in Eq. ()A.20p turn out to be equal to zero 
when this irrep is not contained in the decomposition of the full natural irrep F into its irreducible parts 
Tj. This gives rise to additional selection rules which reduce the number of possible vibrational modes in 
the considered bush. Trying every irrep Tj in Eq. ()A.2ip and analyzing the above mentioned decomposition 
of the natural representation F, we obtain explicitly the whole bush of modes with the symmetry group G. 
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